diff --git a/Cargo.toml b/Cargo.toml index a1044c0..9ba3997 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "fixed_analytics" -version = "2.1.0" +version = "3.0.0" edition = "2024" rust-version = "1.95" authors = ["David Gathercole"] @@ -23,7 +23,7 @@ verify-no-panic = ["dep:no-panic"] [dependencies] fixed = "1.31" -no-panic = { version = "0.1", optional = true } +no-panic = { version = "0.1.37", optional = true } [dev-dependencies] criterion = { version = "0.8", features = ["html_reports"] } diff --git a/README.md b/README.md index 62f1357..8189e81 100644 --- a/README.md +++ b/README.md @@ -30,14 +30,14 @@ Requires Rust 1.95 or later. ```toml [dependencies] -fixed_analytics = "2.1.0" +fixed_analytics = "3.0.0" ``` For `no_std` environments: ```toml [dependencies] -fixed_analytics = { version = "2.0.1", default-features = false } +fixed_analytics = { version = "3.0.0", default-features = false } ``` ## Available Functions @@ -51,7 +51,7 @@ fixed_analytics = { version = "2.0.1", default-features = false } |----------|-----------------|-------------------| | Trigonometric | `sin`, `cos`, `tan`, `sin_cos`, `atan`, `atan2` | `asin`, `acos` | | Hyperbolic | `sinh`, `cosh`, `tanh`, `sinh_cosh`, `asinh` | `acosh`, `atanh`, `acoth`, `coth` | -| Exponential | `exp`, `pow2` | `ln`, `log2`, `log10` | +| Exponential | `exp`, `pow2` | `ln`, `log2`, `log10`, `pow` | | Algebraic | — | `sqrt` | Functions are calculated via polynomial evaluation, CORDIC, and Newton-Raphson techniques. Complete absence of panic is verified at the linker level via the [`no-panic`](https://github.com/dtolnay/no-panic) crate. @@ -82,21 +82,21 @@ Relative error statistics measured against MPFR reference implementations. Accur | sin | 6.06e-4 | 8.78e-5 | 1.28e-3 | 1.16e-8 | 1.68e-9 | 2.43e-8 | | cos | 6.45e-4 | 9.03e-5 | 1.38e-3 | 1.22e-8 | 1.72e-9 | 2.64e-8 | | tan | 7.20e-5 | 3.57e-5 | 2.20e-4 | 1.28e-9 | 3.98e-10 | 3.03e-9 | -| asin | 1.85e-4 | 3.60e-5 | 4.75e-4 | 4.49e-9 | 7.35e-10 | 8.89e-9 | -| acos | 2.33e-5 | 1.47e-5 | 6.86e-5 | 4.51e-10 | 2.71e-10 | 1.38e-9 | -| atan | 1.50e-5 | 1.23e-5 | 3.45e-5 | 3.01e-10 | 2.44e-10 | 7.05e-10 | +| asin | 1.13e-4 | 2.80e-5 | 3.83e-4 | 3.68e-9 | 5.72e-10 | 6.40e-9 | +| acos | 1.79e-5 | 1.17e-5 | 4.81e-5 | 3.50e-10 | 2.11e-10 | 1.04e-9 | +| atan | 7.79e-6 | 6.49e-6 | 1.82e-5 | 1.88e-10 | 1.51e-10 | 4.38e-10 | | sinh | 9.80e-5 | 6.23e-5 | 2.79e-4 | 1.52e-9 | 9.64e-10 | 4.29e-9 | | cosh | 9.40e-5 | 5.75e-5 | 2.77e-4 | 1.44e-9 | 8.90e-10 | 4.25e-9 | | tanh | 1.60e-5 | 1.32e-5 | 2.56e-5 | 2.25e-10 | 1.22e-10 | 3.90e-10 | | coth | 6.68e-6 | 3.54e-6 | 1.80e-5 | 1.41e-10 | 1.16e-10 | 2.74e-10 | -| asinh | 6.44e-4 | 4.83e-4 | 1.75e-3 | 1.03e-8 | 7.59e-9 | 2.85e-8 | -| acosh | 6.74e-4 | 5.21e-4 | 1.80e-3 | 1.05e-8 | 7.96e-9 | 2.88e-8 | -| atanh | 3.01e-4 | 5.90e-5 | 6.25e-4 | 6.68e-9 | 1.32e-9 | 1.44e-8 | -| acoth | 2.10e-3 | 1.33e-3 | 6.67e-3 | 4.26e-8 | 2.62e-8 | 1.39e-7 | +| asinh | 2.42e-5 | 1.61e-5 | 5.08e-5 | 6.27e-10 | 5.01e-10 | 1.19e-9 | +| acosh | 1.87e-5 | 1.48e-5 | 4.65e-5 | 5.29e-10 | 4.75e-10 | 1.11e-9 | +| atanh | 2.21e-4 | 3.99e-5 | 3.33e-4 | 3.63e-9 | 7.20e-10 | 6.24e-9 | +| acoth | 1.21e-3 | 6.85e-4 | 4.09e-3 | 1.94e-8 | 1.23e-8 | 6.28e-8 | | exp | 5.98e-3 | 1.47e-5 | 4.13e-2 | 9.49e-8 | 1.50e-9 | 6.50e-7 | -| ln | 1.35e-5 | 8.76e-6 | 2.97e-5 | 4.50e-10 | 3.48e-10 | 9.17e-10 | -| log2 | 1.33e-5 | 8.48e-6 | 2.92e-5 | 3.46e-10 | 2.24e-10 | 7.21e-10 | -| log10 | 1.44e-5 | 9.28e-6 | 3.14e-5 | 4.49e-10 | 3.27e-10 | 9.07e-10 | +| ln | 1.18e-5 | 7.61e-6 | 2.09e-5 | 4.35e-10 | 3.78e-10 | 6.47e-10 | +| log2 | 9.98e-6 | 6.72e-6 | 2.01e-5 | 1.92e-10 | 1.31e-10 | 3.90e-10 | +| log10 | 1.25e-5 | 8.96e-6 | 2.36e-5 | 4.06e-10 | 3.57e-10 | 6.39e-10 | | pow2 | 3.62e-4 | 2.24e-5 | 2.37e-3 | 5.64e-9 | 4.29e-10 | 3.67e-8 | -| sqrt | 1.77e-7 | 1.16e-7 | 4.74e-7 | 2.70e-12 | 1.78e-12 | 7.16e-12 | +| sqrt | 8.88e-8 | 5.80e-8 | 2.42e-7 | 1.37e-12 | 8.85e-13 | 3.62e-12 | \ No newline at end of file diff --git a/benches/benchmarks.rs b/benches/benchmarks.rs index a8786f5..399fa3a 100644 --- a/benches/benchmarks.rs +++ b/benches/benchmarks.rs @@ -1,68 +1,84 @@ -//! Benchmarks for CORDIC functions. +//! Benchmarks for CORDIC functions, on `I16F16` and on `I64F64`. #![allow(missing_docs, reason = "benchmark code does not need documentation")] use std::hint::black_box; use criterion::{Criterion, criterion_group, criterion_main}; -use fixed::types::I16F16; +use fixed::types::{I16F16, I64F64}; use fixed_analytics::{ - acos, acosh, acoth, asin, asinh, atan, atan2, atanh, cos, cosh, coth, exp, ln, log2, log10, - sin, sin_cos, sinh, sinh_cosh, sqrt, tan, tanh, + CordicNumber, acos, acosh, acoth, asin, asinh, atan, atan2, atanh, cos, cosh, coth, exp, ln, + log2, log10, pow, pow2, sin, sin_cos, sinh, sinh_cosh, sqrt, tan, tanh, }; -fn bench_circular(c: &mut Criterion) { - let angle = I16F16::from_num(0.5); - let x = I16F16::from_num(0.5); +fn bench_type(c: &mut Criterion, name: &str) { + let angle = T::from_num(0.5); + let large_angle = T::from_num(1000.0); + let x = T::from_num(0.5); + let large_x = T::from_num(1.5); + let pos_x = T::from_num(2.0); + let big = T::from_num(50.0); + let small = T::from_num(0.005); - c.bench_function("sin", |b| b.iter(|| sin(black_box(angle)))); - c.bench_function("cos", |b| b.iter(|| cos(black_box(angle)))); - c.bench_function("tan", |b| b.iter(|| tan(black_box(angle)))); - c.bench_function("sin_cos", |b| b.iter(|| sin_cos(black_box(angle)))); - c.bench_function("asin", |b| b.iter(|| asin(black_box(x)))); - c.bench_function("acos", |b| b.iter(|| acos(black_box(x)))); - c.bench_function("atan", |b| b.iter(|| atan(black_box(x)))); - c.bench_function("atan2", |b| { - b.iter(|| atan2(black_box(x), black_box(I16F16::ONE))); - }); + { + let mut g = c.benchmark_group(format!("{name}/circular")); + g.bench_function("sin", |b| b.iter(|| sin(black_box(angle)))); + g.bench_function("cos", |b| b.iter(|| cos(black_box(angle)))); + g.bench_function("tan", |b| b.iter(|| tan(black_box(angle)))); + g.bench_function("sin_cos", |b| b.iter(|| sin_cos(black_box(angle)))); + g.bench_function("sin_cos_large", |b| { + b.iter(|| sin_cos(black_box(large_angle))); + }); + g.bench_function("asin", |b| b.iter(|| asin(black_box(x)))); + g.bench_function("acos", |b| b.iter(|| acos(black_box(x)))); + g.bench_function("atan", |b| b.iter(|| atan(black_box(x)))); + g.bench_function("atan2", |b| { + b.iter(|| atan2(black_box(x), black_box(T::one()))); + }); + g.finish(); + } + { + let mut g = c.benchmark_group(format!("{name}/hyperbolic")); + g.bench_function("sinh", |b| b.iter(|| sinh(black_box(x)))); + g.bench_function("cosh", |b| b.iter(|| cosh(black_box(x)))); + g.bench_function("tanh", |b| b.iter(|| tanh(black_box(x)))); + g.bench_function("coth", |b| b.iter(|| coth(black_box(x)))); + g.bench_function("sinh_cosh", |b| b.iter(|| sinh_cosh(black_box(x)))); + g.bench_function("asinh", |b| b.iter(|| asinh(black_box(x)))); + g.bench_function("asinh_large", |b| b.iter(|| asinh(black_box(big)))); + g.bench_function("acosh", |b| b.iter(|| acosh(black_box(large_x)))); + g.bench_function("atanh", |b| b.iter(|| atanh(black_box(x)))); + g.bench_function("acoth", |b| b.iter(|| acoth(black_box(large_x)))); + g.finish(); + } + { + let mut g = c.benchmark_group(format!("{name}/exponential")); + g.bench_function("exp", |b| b.iter(|| exp(black_box(x)))); + g.bench_function("pow2", |b| b.iter(|| pow2(black_box(x)))); + g.bench_function("ln", |b| b.iter(|| ln(black_box(pos_x)))); + g.bench_function("log2", |b| b.iter(|| log2(black_box(pos_x)))); + g.bench_function("log10", |b| b.iter(|| log10(black_box(pos_x)))); + g.bench_function("pow", |b| { + b.iter(|| pow(black_box(pos_x), black_box(x))); + }); + g.finish(); + } + { + let mut g = c.benchmark_group(format!("{name}/algebraic")); + g.bench_function("sqrt", |b| b.iter(|| sqrt(black_box(pos_x)))); + g.bench_function("sqrt_large", |b| b.iter(|| sqrt(black_box(big)))); + g.bench_function("sqrt_small", |b| b.iter(|| sqrt(black_box(small)))); + g.finish(); + } } -fn bench_hyperbolic(c: &mut Criterion) { - let x = I16F16::from_num(0.5); - let large_x = I16F16::from_num(1.5); - - c.bench_function("sinh", |b| b.iter(|| sinh(black_box(x)))); - c.bench_function("cosh", |b| b.iter(|| cosh(black_box(x)))); - c.bench_function("tanh", |b| b.iter(|| tanh(black_box(x)))); - c.bench_function("coth", |b| b.iter(|| coth(black_box(x)))); - c.bench_function("sinh_cosh", |b| b.iter(|| sinh_cosh(black_box(x)))); - c.bench_function("asinh", |b| b.iter(|| asinh(black_box(x)))); - c.bench_function("acosh", |b| b.iter(|| acosh(black_box(large_x)))); - c.bench_function("atanh", |b| b.iter(|| atanh(black_box(x)))); - c.bench_function("acoth", |b| b.iter(|| acoth(black_box(large_x)))); +fn bench_i16f16(c: &mut Criterion) { + bench_type::(c, "I16F16"); } -fn bench_exponential(c: &mut Criterion) { - let x = I16F16::from_num(0.5); - let pos_x = I16F16::from_num(2.0); - - c.bench_function("exp", |b| b.iter(|| exp(black_box(x)))); - c.bench_function("ln", |b| b.iter(|| ln(black_box(pos_x)))); - c.bench_function("log2", |b| b.iter(|| log2(black_box(pos_x)))); - c.bench_function("log10", |b| b.iter(|| log10(black_box(pos_x)))); -} - -fn bench_algebraic(c: &mut Criterion) { - let x = I16F16::from_num(2.0); - - c.bench_function("sqrt", |b| b.iter(|| sqrt(black_box(x)))); +fn bench_i64f64(c: &mut Criterion) { + bench_type::(c, "I64F64"); } -criterion_group!( - benches, - bench_circular, - bench_hyperbolic, - bench_exponential, - bench_algebraic -); +criterion_group!(benches, bench_i16f16, bench_i64f64); criterion_main!(benches); diff --git a/codecov.yml b/codecov.yml index 527f36d..521e421 100644 --- a/codecov.yml +++ b/codecov.yml @@ -20,4 +20,5 @@ comment: behavior: default ignore: - - "benches/**/*" + - "benches" + - "src/bin" diff --git a/src/bin/verify_no_panic.rs b/src/bin/verify_no_panic.rs index 6d1c735..47d38c9 100644 --- a/src/bin/verify_no_panic.rs +++ b/src/bin/verify_no_panic.rs @@ -1,26 +1,25 @@ -//! Binary that instantiates every public function with a concrete type. +//! Binary that instantiates every public function with concrete types. //! //! This exists solely to trigger monomorphization so that `no_panic`'s //! linker-level check can verify that no panic paths survive optimization. -//! It is only compiled under the `verify-no-panic` feature. +//! It is only compiled under the `verify-no-panic` feature, for three +//! layouts that compile to different code: `I16F16`, `I24F8` (more integer +//! than fractional bits) and `I64F64` (128-bit). #[cfg(not(feature = "verify-no-panic"))] compile_error!("this binary should only be built with --features verify-no-panic"); -use fixed::types::I16F16; +use fixed::types::{I16F16, I24F8, I64F64}; +use fixed_analytics::CordicNumber; use fixed_analytics::bounded::{NonNegative, OpenUnitInterval}; use fixed_analytics::ops::algebraic::sqrt_nonneg; use fixed_analytics::ops::hyperbolic::atanh_open; use fixed_analytics::{ acos, acosh, acoth, asin, asinh, atan, atan2, atanh, cos, cosh, coth, exp, ln, log2, log10, - pow2, sin, sin_cos, sinh, sinh_cosh, sqrt, tan, tanh, + pow, pow2, sin, sin_cos, sinh, sinh_cosh, sqrt, tan, tanh, }; -fn main() { - // Use black_box to prevent the optimizer from eliminating calls entirely. - let x = std::hint::black_box(I16F16::from_num(0.5)); - let y = std::hint::black_box(I16F16::from_num(0.25)); - +fn exercise(x: T, y: T, two: T) { // Total functions (return T) let _ = std::hint::black_box(sin(x)); let _ = std::hint::black_box(cos(x)); @@ -35,6 +34,7 @@ fn main() { let _ = std::hint::black_box(tanh(x)); let _ = std::hint::black_box(sinh_cosh(x)); let _ = std::hint::black_box(asinh(x)); + let _ = std::hint::black_box(asinh(two)); // Fallible functions (return Result) let _ = std::hint::black_box(asin(x)); @@ -43,14 +43,36 @@ fn main() { let _ = std::hint::black_box(ln(x)); let _ = std::hint::black_box(log2(x)); let _ = std::hint::black_box(log10(x)); - let _ = std::hint::black_box(acosh(I16F16::from_num(2))); + let _ = std::hint::black_box(acosh(two)); let _ = std::hint::black_box(atanh(x)); let _ = std::hint::black_box(coth(x)); - let _ = std::hint::black_box(acoth(I16F16::from_num(2))); + let _ = std::hint::black_box(acoth(two)); + let _ = std::hint::black_box(pow(two, x)); // Type-safe wrapper functions - let nn = NonNegative::new(x).unwrap(); - let _ = std::hint::black_box(sqrt_nonneg(nn)); - let ou = OpenUnitInterval::new(x).unwrap(); - let _ = std::hint::black_box(atanh_open(ou)); + if let Some(nn) = NonNegative::new(x) { + let _ = std::hint::black_box(sqrt_nonneg(nn)); + } + if let Some(ou) = OpenUnitInterval::new(x) { + let _ = std::hint::black_box(atanh_open(ou)); + } +} + +fn main() { + // Use black_box to prevent the optimizer from eliminating calls entirely. + exercise( + std::hint::black_box(I16F16::from_num(0.5)), + std::hint::black_box(I16F16::from_num(0.25)), + std::hint::black_box(I16F16::from_num(2)), + ); + exercise( + std::hint::black_box(I24F8::from_num(0.5)), + std::hint::black_box(I24F8::from_num(0.25)), + std::hint::black_box(I24F8::from_num(2)), + ); + exercise( + std::hint::black_box(I64F64::from_num(0.5)), + std::hint::black_box(I64F64::from_num(0.25)), + std::hint::black_box(I64F64::from_num(2)), + ); } diff --git a/src/bounded.rs b/src/bounded.rs index a833970..7da2ecd 100644 --- a/src/bounded.rs +++ b/src/bounded.rs @@ -103,6 +103,13 @@ impl UnitInterval { (value >= -one && value <= one).then_some(Self(value)) } + /// Constructs `1 / x` where `x >= 1`, which is always in (0, 1]. + #[inline] + #[must_use] + pub fn from_reciprocal(x: AtLeastOne) -> Self { + Self(T::one().div(x.0)) + } + /// Unwraps the inner value. #[inline] #[must_use] @@ -259,6 +266,15 @@ mod tests { assert!(UnitInterval::new(I16F16::from_num(-1.1)).is_none()); } + #[test] + fn unit_interval_from_reciprocal() { + let at_least = AtLeastOne::new(I16F16::from_num(4)).unwrap(); + assert_eq!( + UnitInterval::from_reciprocal(at_least).get(), + I16F16::from_num(0.25) + ); + } + #[test] fn unit_interval_get() { let unit = UnitInterval::new(I16F16::from_num(0.5)).unwrap(); diff --git a/src/kernel/cordic.rs b/src/kernel/cordic.rs index b1f42a9..f193c6e 100644 --- a/src/kernel/cordic.rs +++ b/src/kernel/cordic.rs @@ -20,25 +20,28 @@ //! - σ = ±1 (direction of rotation) //! - d = +1 for circular, -1 for hyperbolic, 0 for linear //! - angle[i] = atan(2^-i) for circular, atanh(2^-i) for hyperbolic +//! +//! Both kernels stop after about a third of the fractional bits and close +//! the residual angle with one division: for `|w| = |y/x| < 2^-⌈frac/3⌉`, +//! `atan(w)` and `atanh(w)` differ from `w` by under `|w|³/3 < 2^-frac/3`. use crate::tables::hyperbolic::needs_repeat; use crate::tables::{ATAN_TABLE, ATANH_TABLE}; use crate::traits::CordicNumber; -/// Table lookup for CORDIC iteration. -/// -/// Index is bounded by CORDIC iteration limits: -/// - Circular mode: `min(frac_bits, 62)` → max index 61 -/// - Hyperbolic mode: `min(frac_bits, 54)` with `i.saturating_sub(1)` → max index 53 -/// -/// Since the tables have 64 elements and max index is 61, bounds are always satisfied. +/// Index of the last shift stage before the closing division: after the +/// stage with shift `2^-k` the residual angle, and so `|y/x|`, is below +/// `2^-k`. At most 42, so every table lookup is in bounds. +const fn last_stage(frac_bits: u32) -> u32 { + frac_bits.div_ceil(3) +} + +/// Table lookup for CORDIC iteration. The modulo is a no-op that keeps the +/// access provably in bounds for the no-panic check. #[inline] const fn table_lookup(table: &[i64; 64], index: u32) -> i64 { - #[allow( - clippy::indexing_slicing, - reason = "index bounded by CORDIC iteration limits" - )] - table[index as usize] + #[allow(clippy::indexing_slicing, reason = "index reduced modulo table length")] + table[index as usize % table.len()] } /// Converts an I1F63 table constant to `T`, rounding to nearest. @@ -73,21 +76,21 @@ fn from_i1f63_rounded(bits: i64) -> T { /// Performs circular CORDIC in vectoring mode. /// -/// Given an initial vector (x, y), rotates it until y ≈ 0. -/// After iteration: -/// - x ≈ K * sqrt(x₀² + y₀²) -/// - y ≈ 0 -/// - z ≈ z₀ + atan(y₀/x₀) +/// Given an initial vector (x, y) with x > 0, rotates it toward the x axis +/// through `⌈frac/3⌉ + 1` shift stages, then closes the residual angle +/// with one division. /// /// # Arguments /// -/// * `x` - Initial x coordinate (should be positive for standard use) +/// * `x` - Initial x coordinate (must be positive) /// * `y` - Initial y coordinate /// * `z` - Initial angle accumulator (usually 0) /// /// # Returns /// -/// Tuple of (x, y, z) after CORDIC iterations. +/// Tuple of (x, y, z): z ≈ z₀ + atan(y₀/x₀); x ≈ K·sqrt(x₀² + y₀²) with K +/// the gain of the stages run; y is the residual the shift stages left, +/// whose angle is already in z. /// /// # Note /// @@ -95,9 +98,8 @@ fn from_i1f63_rounded(bits: i64) -> T { #[must_use] pub fn circular_vectoring(mut x: T, mut y: T, mut z: T) -> (T, T, T) { let zero = T::zero(); - let iterations = T::frac_bits().min(62); - for i in 0..iterations { + for i in 0..=last_stage(T::frac_bits()) { let angle = from_i1f63_rounded::(table_lookup(&ATAN_TABLE, i)); if y < zero { @@ -115,17 +117,17 @@ pub fn circular_vectoring(mut x: T, mut y: T, mut z: T) -> (T, } } + // Close the residual angle: atan(y/x) = y/x to within |y/x|³/3. + z = z.saturating_add(y.div(x)); + (x, y, z) } /// Performs hyperbolic CORDIC in vectoring mode. /// -/// Drives y toward zero while accumulating the hyperbolic angle. -/// -/// After iteration: -/// - x ≈ `K_h` * sqrt(x₀² - y₀²) (for |x| > |y|) -/// - y ≈ 0 -/// - z ≈ z₀ + atanh(y₀/x₀) +/// Drives y toward zero through shift stages 1..=⌈frac/3⌉ (repeating 4, +/// 13, 40, … for convergence), then closes the residual angle with one +/// division. /// /// # Arguments /// @@ -135,7 +137,9 @@ pub fn circular_vectoring(mut x: T, mut y: T, mut z: T) -> (T, /// /// # Returns /// -/// Tuple of (x, y, z) after CORDIC iterations. +/// Tuple of (x, y, z): z ≈ z₀ + atanh(y₀/x₀); x ≈ `K_h`·sqrt(x₀² - y₀²) +/// with `K_h` the gain of the stages run; y is the residual the shift +/// stages left, whose angle is already in z. /// /// # Note /// @@ -144,52 +148,36 @@ pub fn circular_vectoring(mut x: T, mut y: T, mut z: T) -> (T, #[must_use] pub fn hyperbolic_vectoring(mut x: T, mut y: T, mut z: T) -> (T, T, T) { let zero = T::zero(); - // Use at least 24 iterations for better accuracy, even for lower precision types. - let max_iterations = T::frac_bits().clamp(24, 54); - - let mut i: u32 = 1; - let mut iteration_count: u32 = 0; - let mut repeated = false; - - while iteration_count < max_iterations && i < 64 { - let table_index = i.saturating_sub(1); - let angle = T::from_i1f63(table_lookup(&ATANH_TABLE, table_index)); - - // Hyperbolic pseudo-rotation equations: - // x' = x + σ*y*2^(-i) - // y' = y + σ*x*2^(-i) - // z' = z + σ*angle (accumulating for vectoring) - // where σ = -sign(y) to drive y toward zero - - if y < zero { - // y is negative: σ = +1 - // x' = x + y*2^(-i) [y is negative, so this subtracts magnitude] - // y' = y + x*2^(-i) [adds positive to make less negative] - // z' = z - angle [accumulate negative contribution] - let x_new = x.saturating_add(y >> i); - y = y.saturating_add(x >> i); - x = x_new; - z -= angle; - } else { - // y is positive or zero: σ = -1 - // x' = x - y*2^(-i) [subtracts positive] - // y' = y - x*2^(-i) [subtracts to decrease toward zero] - // z' = z + angle [accumulate positive contribution] - let x_new = x.saturating_sub(y >> i); - y = y.saturating_sub(x >> i); - x = x_new; - z += angle; - } - iteration_count += 1; + for i in 1..=last_stage(T::frac_bits()) { + let angle = T::from_i1f63(table_lookup(&ATANH_TABLE, i - 1)); + // Stages 4, 13, 40, … run twice so the sequence converges. + let passes = if needs_repeat(i) { 2 } else { 1 }; - if needs_repeat(i) && !repeated { - repeated = true; - } else { - repeated = false; - i += 1; + for _ in 0..passes { + // Hyperbolic pseudo-rotation equations: + // x' = x + σ*y*2^(-i) + // y' = y + σ*x*2^(-i) + // z' = z + σ*angle (accumulating for vectoring) + // where σ = -sign(y) to drive y toward zero + if y < zero { + // y is negative: σ = +1 + let x_new = x.saturating_add(y >> i); + y = y.saturating_add(x >> i); + x = x_new; + z -= angle; + } else { + // y is positive or zero: σ = -1 + let x_new = x.saturating_sub(y >> i); + y = y.saturating_sub(x >> i); + x = x_new; + z += angle; + } } } + // Close the residual angle: atanh(y/x) = y/x to within |y/x|³/3. + z = z.saturating_add(y.div(x)); + (x, y, z) } diff --git a/src/lib.rs b/src/lib.rs index 18d3d45..f667fd6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -27,7 +27,7 @@ //! |--------------|-------|----------| //! | Trigonometric | [`sin`], [`cos`], [`tan`], [`sin_cos`], [`atan`], [`atan2`] | [`asin`], [`acos`] | //! | Hyperbolic | [`sinh`], [`cosh`], [`tanh`], [`sinh_cosh`], [`asinh`] | [`acosh`], [`atanh`], [`acoth`], [`coth`] | -//! | Exponential | [`exp`], [`pow2`] | [`ln`], [`log2`], [`log10`] | +//! | Exponential | [`exp`], [`pow2`] | [`ln`], [`log2`], [`log10`], [`pow`] | //! | Algebraic | — | [`sqrt`] | //! //! Functions use polynomial evaluation, CORDIC, and Newton-Raphson techniques. @@ -70,5 +70,5 @@ pub use traits::CordicNumber; // Re-export all mathematical functions at crate root for convenience pub use ops::algebraic::sqrt; pub use ops::circular::{acos, asin, atan, atan2, cos, sin, sin_cos, tan}; -pub use ops::exponential::{exp, ln, log2, log10, pow2}; +pub use ops::exponential::{exp, ln, log2, log10, pow, pow2}; pub use ops::hyperbolic::{acosh, acoth, asinh, atanh, cosh, coth, sinh, sinh_cosh, tanh}; diff --git a/src/ops/algebraic.rs b/src/ops/algebraic.rs index 2647aeb..5c17890 100644 --- a/src/ops/algebraic.rs +++ b/src/ops/algebraic.rs @@ -4,7 +4,7 @@ use crate::bounded::NonNegative; use crate::error::{Error, Result}; use crate::traits::CordicNumber; -/// Square root. Domain: `x ≥ 0`. Uses Newton-Raphson iteration. +/// Square root, rounded to the nearest representable value. Domain: `x ≥ 0`. /// /// # Errors /// Returns `DomainError` if `x < 0`. @@ -16,91 +16,20 @@ pub fn sqrt(x: T) -> Result { .ok_or_else(|| Error::domain("sqrt", "non-negative value")) } -/// Infallible square root for non-negative values. +/// Infallible square root for non-negative values, rounded to the nearest +/// representable value. /// /// This function takes a [`NonNegative`] wrapper, guaranteeing at the type /// level that the input is valid. No domain check is performed at runtime. /// /// Use this when the non-negativity of the input is already established /// through mathematical invariants (e.g., `1 + x²`, `1 - x²` for `|x| ≤ 1`). +/// +/// The result is correctly rounded at every magnitude (raw bits +/// `⌊√(X·2^f) + ½⌋`), computed from an integer square root of the raw +/// representation; see [`CordicNumber::sqrt_round`]. #[must_use] #[cfg_attr(feature = "verify-no-panic", no_panic::no_panic)] pub fn sqrt_nonneg(x: NonNegative) -> T { - let x = x.get(); - let zero = T::zero(); - let one = T::one(); - let half = T::half(); - - if x == zero { - return zero; - } - - if x == one { - return one; - } - - // Initial guess: use bit-level estimation for faster convergence - // For sqrt(x), a good initial estimate is 2^(floor(log2(x))/2) - let mut guess = if x > one { - // For large numbers, use bit estimation: sqrt(x) ≈ 2^(log2(x)/2) - // We find the approximate position by successive squaring comparison - let mut g = one; - let four = T::two().saturating_mul(T::two()); - let mut test = x; - // Guard: max 64 iterations sufficient for any representable value - let mut iter_guard = 0u32; - while test >= four && iter_guard < 64 { - test = test >> 2; - g = g << 1; - iter_guard += 1; - } - // g is now approximately sqrt(x) rounded down to a power of 2 - // Refine: average with x/g for a better starting point - let quotient = x.div(g); - g.saturating_add(quotient) >> 1 - } else { - // For numbers < 1, x is a reasonable starting point - // since sqrt(x) > x when 0 < x < 1 - x - }; - - // Newton-Raphson iteration: x_new = (x_old + n/x_old) / 2 - // Number of iterations depends on precision needed. - // Newton-Raphson for sqrt converges quadratically, so 8-20 iterations - // is sufficient for any fixed-point precision up to 128 bits. - let iterations = (T::frac_bits() / 2).clamp(8, 20); - - // Pre-compute epsilon: approximately 2^(-frac_bits/2) for convergence check. - // frac_bits ≤ 128 for all supported types, so shift is in range [0, 63]. - #[allow( - clippy::cast_possible_wrap, - clippy::cast_sign_loss, - reason = "frac_bits bounded by type size" - )] - let epsilon_shift = (63i32 - (T::frac_bits() / 2) as i32).max(0) as u32; - let epsilon = T::from_i1f63(1i64 << epsilon_shift); - - // Run iterations - 1 times with early exit on convergence - for _ in 0..iterations.saturating_sub(1) { - let quotient = x.div(guess); - let sum = guess.saturating_add(quotient); - let new_guess = sum.saturating_mul(half); - - let diff = if new_guess > guess { - new_guess.saturating_sub(guess) - } else { - guess.saturating_sub(new_guess) - }; - - if diff <= epsilon { - return new_guess; - } - - guess = new_guess; - } - - // Final iteration - always performed, result always returned - let quotient = x.div(guess); - let sum = guess.saturating_add(quotient); - sum.saturating_mul(half) + x.get().sqrt_round() } diff --git a/src/ops/circular.rs b/src/ops/circular.rs index 8e810a4..153939c 100644 --- a/src/ops/circular.rs +++ b/src/ops/circular.rs @@ -7,6 +7,9 @@ use crate::ops::algebraic::sqrt_nonneg; use crate::tables::chebyshev::{COS_Q_HI, COS_Q_LO, SIN_P_HI, SIN_P_LO, horner}; use crate::traits::CordicNumber; +/// `1/(2π)` as I1F63, for the angle reduction multiply in [`sin_cos`]. +const FRAC_1_2PI_I1F63: i64 = 0x145F_306D_C9C8_82A5; + /// Sine and cosine. More efficient than separate calls. Accepts any angle. #[must_use] #[cfg_attr(feature = "verify-no-panic", no_panic::no_panic)] @@ -18,10 +21,20 @@ pub fn sin_cos(angle: T) -> (T, T) { // Reduce angle to [-π, π] using direct quotient computation. // This handles arbitrarily large angles without iteration limits. let reduced = if angle > pi || angle < -pi { - // Compute n = round(angle / 2π), then reduced = angle - n * 2π - let quotient = angle.div(two_pi); + // Compute n = round(angle / 2π), then reduced = angle - n * 2π. + // The reciprocal multiply is within |angle|·1.5·2^-frac + 2^-frac < 1 + // of the true quotient when the integer bits do not outnumber the + // fractional bits, so n is off by at most one and the clamp below + // absorbs it bit for bit; other types keep the division. + let quotient = if T::total_bits() <= 2 * T::frac_bits() { + angle.saturating_mul(T::from_i1f63(FRAC_1_2PI_I1F63)) + } else { + angle.div(two_pi) + }; let n = quotient.round(); - angle.saturating_sub(n.saturating_mul(two_pi)) + // n·2π is exact and the difference small, so wrapping arithmetic + // recovers it even when n·2π exceeds MAX (angles within π of MAX). + angle.wrapping_sub(n.wrapping_mul(two_pi)) } else { angle }; @@ -164,6 +177,11 @@ pub fn tan(angle: T) -> T { /// /// # Errors /// Returns `DomainError` if `|x| > 1`. +#[allow( + clippy::inline_always, + reason = "inlined into acos so the no-panic link check sees only nounwind callees" +)] +#[inline(always)] #[must_use = "returns the arcsine result which should be handled"] #[cfg_attr(feature = "verify-no-panic", no_panic::no_panic)] pub fn asin(x: T) -> Result { diff --git a/src/ops/exponential.rs b/src/ops/exponential.rs index c5c30f1..f90ff59 100644 --- a/src/ops/exponential.rs +++ b/src/ops/exponential.rs @@ -5,6 +5,9 @@ use crate::error::{Error, Result}; use crate::ops::hyperbolic::atanh_open; use crate::traits::CordicNumber; +/// `1/ln 2 − 1` as I1F63 (`1/ln 2` itself exceeds 1). +const FRAC_1_LN2_MINUS_1_I1F63: i64 = 0x38AA_3B29_5C17_F0BC; + /// Right-shifts a non-negative value by `n`, rounding to nearest. /// /// A plain `>>` truncates, which biases results low by up to a full ulp. @@ -70,57 +73,79 @@ pub fn exp(x: T) -> T { return one; } - // Argument reduction: exp(x) = 2^k * exp(r), where r ∈ [0, ln2). - // Compute k = floor(x / ln2): start from the truncated quotient and - // adjust once if the remainder is negative (the quotient error is - // below one ulp, so a single adjustment suffices). Keeping r - // non-negative gives exp(r) ∈ [1, 2) — a full significand ahead of - // the final scaling shift. - #[allow(clippy::cast_possible_wrap, reason = "total_bits bounded by type size")] - let max_shift = (T::total_bits() - 1) as i32; - let mut scale = x.div(ln2).to_i32(); + // Argument reduction: exp(x) = 2^k * exp(r), where r ∈ [0, ln2), so + // exp(r) ∈ [1, 2) carries a full significand into the scaling shift. + // k = ⌊x / ln2⌋ is estimated to within one (`to_i32` floors; the + // reciprocal multiply is within |x|·1.5·2^-frac + 2^-frac < 1 when the + // integer bits do not outnumber the fractional bits, and other types + // divide) and then corrected once in either direction. + #[allow(clippy::cast_possible_wrap, reason = "bit counts bounded by type size")] + let (int_bits, frac_bits) = ( + (T::total_bits() - T::frac_bits()) as i32, + T::frac_bits() as i32, + ); + let quotient = if T::total_bits() <= 2 * T::frac_bits() { + x.saturating_add(x.saturating_mul(T::from_i1f63(FRAC_1_LN2_MINUS_1_I1F63))) + } else { + x.div(ln2) + }; + let mut scale = quotient.to_i32(); - // Early exit for values that will saturate after scaling. This also - // keeps `scale` small enough that scale * ln2 below cannot overflow. - if scale > max_shift { + // 2^k·exp(r) exceeds MAX once k ≥ int_bits − 1 and rounds to zero once + // k ≤ −(frac_bits + 2); allow one of slack until k is corrected. + if scale >= int_bits { return T::max_value(); } - if scale < -max_shift { + if scale < -(frac_bits + 2) { return zero; } - let mut r = x.saturating_sub(T::from_num(scale).saturating_mul(ln2)); + let mut r = x.saturating_sub(ln2.mul_int(scale)); if r < zero { scale -= 1; r = r.saturating_add(ln2); - if scale < -max_shift { - return zero; - } + } else if r >= ln2 { + scale += 1; + r = r.saturating_sub(ln2); + } + if scale >= int_bits - 1 { + return T::max_value(); + } + if scale <= -(frac_bits + 2) { + return zero; } - // Factored Taylor: exp(r) = 1 + r*(1 + r/2*(1 + r/3*(1 + ... r/n))) + // Factored Taylor: exp(r) = 1 + r*(1 + r/2*(1 + r/3*(1 + ... r/n))). + // The omitted term (ln2)^(n+1)/(n+1)! is 1.3e-6 at degree 7, 1.4e-12 at + // degree 12 and 2.7e-22 at degree 19: below half an ulp up to 18, 38 + // and 70 fractional bits respectively. let mut p = one; + if T::frac_bits() >= 40 { + p = one.saturating_add(r.div_int(19).saturating_mul(p)); + p = one.saturating_add(r.div_int(18).saturating_mul(p)); + p = one.saturating_add(r.div_int(17).saturating_mul(p)); + p = one.saturating_add(r.div_int(16).saturating_mul(p)); + p = one.saturating_add(r.div_int(15).saturating_mul(p)); + p = one.saturating_add(r.div_int(14).saturating_mul(p)); + p = one.saturating_add(r.div_int(13).saturating_mul(p)); + } if T::frac_bits() >= 24 { - // High precision: degree 12 Taylor - // Truncation error: |r^13/13!| ≤ (ln2)^13/13! ≈ 3.4e-15 - p = one.saturating_add(r.div(T::from_num(12)).saturating_mul(p)); - p = one.saturating_add(r.div(T::from_num(11)).saturating_mul(p)); - p = one.saturating_add(r.div(T::from_num(10)).saturating_mul(p)); - p = one.saturating_add(r.div(T::from_num(9)).saturating_mul(p)); - p = one.saturating_add(r.div(T::from_num(8)).saturating_mul(p)); - } - // Common terms (degree 7 base) - // Low-precision truncation error: |r^8/8!| ≤ (ln2)^8/8! ≈ 8.9e-7 - p = one.saturating_add(r.div(T::from_num(7)).saturating_mul(p)); - p = one.saturating_add(r.div(T::from_num(6)).saturating_mul(p)); - p = one.saturating_add(r.div(T::from_num(5)).saturating_mul(p)); - p = one.saturating_add(r.div(T::from_num(4)).saturating_mul(p)); - p = one.saturating_add(r.div(T::from_num(3)).saturating_mul(p)); - p = one.saturating_add(r.div(T::from_num(2)).saturating_mul(p)); + p = one.saturating_add(r.div_int(12).saturating_mul(p)); + p = one.saturating_add(r.div_int(11).saturating_mul(p)); + p = one.saturating_add(r.div_int(10).saturating_mul(p)); + p = one.saturating_add(r.div_int(9).saturating_mul(p)); + p = one.saturating_add(r.div_int(8).saturating_mul(p)); + } + p = one.saturating_add(r.div_int(7).saturating_mul(p)); + p = one.saturating_add(r.div_int(6).saturating_mul(p)); + p = one.saturating_add(r.div_int(5).saturating_mul(p)); + p = one.saturating_add(r.div_int(4).saturating_mul(p)); + p = one.saturating_add(r.div_int(3).saturating_mul(p)); + p = one.saturating_add(r.div_int(2).saturating_mul(p)); let exp_r = one.saturating_add(r.saturating_mul(p)); // Scale by 2^scale using bit shifts. - // scale is already bounded to [-max_shift, max_shift] by the early exits above. + // scale is already bounded to (-(frac_bits + 2), int_bits - 1) by the exits above. #[allow( clippy::cast_sign_loss, reason = "sign of scale checked before each cast" @@ -141,6 +166,33 @@ pub fn exp(x: T) -> T { } } +/// Power function `base^exponent`, computed as `exp(exponent · ln(base))`. +/// Domain: `base > 0`, or `base = 0` with `exponent ≥ 0`. +/// +/// The logarithm's rounding error is multiplied by the exponent, so the +/// relative error grows with `|exponent|`. Saturates like [`exp`]. +/// +/// # Errors +/// Returns `DomainError` if `base < 0`, or if `base = 0` and `exponent < 0`. +#[must_use = "returns the power result which should be handled"] +#[cfg_attr(feature = "verify-no-panic", no_panic::no_panic)] +pub fn pow(base: T, exponent: T) -> Result { + let zero = T::zero(); + if base < zero || (base == zero && exponent < zero) { + return Err(Error::domain( + "pow", + "positive base, or zero base with non-negative exponent", + )); + } + if exponent == zero { + return Ok(T::one()); + } + if base == zero { + return Ok(zero); + } + Ok(exp(exponent.saturating_mul(ln_positive(base)))) +} + /// Natural logarithm. Domain: `x > 0`. /// /// # Errors @@ -148,16 +200,20 @@ pub fn exp(x: T) -> T { #[must_use = "returns the natural logarithm result which should be handled"] #[cfg_attr(feature = "verify-no-panic", no_panic::no_panic)] pub fn ln(x: T) -> Result { + if x <= T::zero() { + return Err(Error::domain("ln", "positive value")); + } + Ok(ln_positive(x)) +} + +/// Natural logarithm of a positive value. The caller must ensure `x > 0`. +pub(crate) fn ln_positive(x: T) -> T { let zero = T::zero(); let one = T::one(); let two = T::two(); - if x <= zero { - return Err(Error::domain("ln", "positive value")); - } - if x == one { - return Ok(zero); + return zero; } // For x far from 1, use argument reduction: @@ -165,27 +221,28 @@ pub fn ln(x: T) -> Result { // where k is chosen so that x * 2^(-k) is close to 1 let ln2 = T::ln_2(); - let mut normalized = x; - let mut k_ln2 = zero; - - // Reduce to range [0.5, 2] for better convergence. - let half = T::half(); - // For large x, divide by 2 repeatedly - let mut i = 0; - while normalized > two && i < 128 { - normalized = normalized >> 1; - k_ln2 = k_ln2.saturating_add(ln2); - i += 1; - } - - // For small x (< 0.5), multiply by 2 repeatedly - i = 0; - while normalized < half && i < 128 { - normalized = normalized.saturating_add(normalized); - k_ln2 = k_ln2.saturating_sub(ln2); - i += 1; - } + // Reduce to [0.5, 2] from the leading bit e = ⌊log₂ x⌋: x > 2 shifts + // right until ≤ 2 (stopping early when a shift lands exactly on 2), + // x < 0.5 doubles until ≥ 0.5. + let e = x.checked_int_log2().unwrap_or(0); + #[allow( + clippy::cast_sign_loss, + reason = "shift counts are non-negative by construction" + )] + let (normalized, k) = if e >= 1 { + let candidate = x >> (e - 1) as u32; + if candidate == two { + (candidate, e - 1) + } else { + (x >> e as u32, e) + } + } else if e <= -2 { + (x << (-1 - e) as u32, e + 1) + } else { + (x, 0) + }; + let k_ln2 = ln2.mul_int(k); // Now compute ln(normalized) where 0.5 <= normalized <= 2 // Using ln(x) = 2 * atanh((x-1)/(x+1)) @@ -199,7 +256,7 @@ pub fn ln(x: T) -> Result { let atanh_val = atanh_open(arg); let ln_normalized = atanh_val.saturating_add(atanh_val); // 2 * atanh - Ok(ln_normalized.saturating_add(k_ln2)) + ln_normalized.saturating_add(k_ln2) } /// Base-2 logarithm. Domain: `x > 0`. diff --git a/src/ops/hyperbolic.rs b/src/ops/hyperbolic.rs index fe175e1..54df348 100644 --- a/src/ops/hyperbolic.rs +++ b/src/ops/hyperbolic.rs @@ -1,9 +1,10 @@ //! Hyperbolic functions via hyperbolic CORDIC. -use crate::bounded::{AtLeastOne, NonNegative, OpenUnitInterval}; +use crate::bounded::{AtLeastOne, NonNegative, OpenUnitInterval, UnitInterval}; use crate::error::{Error, Result}; use crate::kernel::hyperbolic_vectoring; use crate::ops::algebraic::sqrt_nonneg; +use crate::ops::exponential::ln_positive; use crate::traits::CordicNumber; /// Hyperbolic CORDIC converges for |x| < sum of atanh table ≈ 1.1182. @@ -64,44 +65,38 @@ pub fn sinh_cosh(x: T) -> (T, T) { // step's multiplicand u/K below 1 for K ≥ 2. let u = reduced.saturating_mul(reduced); - // High-precision path requires enough integer bits for divisors up to 182. + // Divisors are (2k)(2k+1) for sinh and (2k-1)(2k) for cosh. At |x| ≤ + // 1.1182 the omitted relative term is 7.7e-8 / 8.0e-9 at degrees 9 / 10, + // 3.7e-12 / 2.9e-13 at 13 / 14 and 4.5e-22 / 2.4e-23 at 21 / 22: below + // half an ulp up to 22, 37 and 70 fractional bits respectively. let mut sp = one; - let (mut sh, mut ch) = if T::frac_bits() >= 24 && T::total_bits() >= T::frac_bits() + 9 { - // High precision: degree 13 sinh, degree 14 cosh - sp = one.saturating_add(u.div(T::from_num(156)).saturating_mul(sp)); - sp = one.saturating_add(u.div(T::from_num(110)).saturating_mul(sp)); - sp = one.saturating_add(u.div(T::from_num(72)).saturating_mul(sp)); - sp = one.saturating_add(u.div(T::from_num(42)).saturating_mul(sp)); - sp = one.saturating_add(u.div(T::from_num(20)).saturating_mul(sp)); - sp = one.saturating_add(u.div(T::from_num(6)).saturating_mul(sp)); - let sinh_approx = reduced.saturating_mul(sp); - - let mut cp = one; - cp = one.saturating_add(u.div(T::from_num(182)).saturating_mul(cp)); - cp = one.saturating_add(u.div(T::from_num(132)).saturating_mul(cp)); - cp = one.saturating_add(u.div(T::from_num(90)).saturating_mul(cp)); - cp = one.saturating_add(u.div(T::from_num(56)).saturating_mul(cp)); - cp = one.saturating_add(u.div(T::from_num(30)).saturating_mul(cp)); - cp = one.saturating_add(u.div(T::from_num(12)).saturating_mul(cp)); - cp = one.saturating_add(u.div(T::from_num(2)).saturating_mul(cp)); - (sinh_approx, cp) - } else { - // Low precision: degree 9 sinh, degree 10 cosh - sp = one; - sp = one.saturating_add(u.div(T::from_num(72)).saturating_mul(sp)); - sp = one.saturating_add(u.div(T::from_num(42)).saturating_mul(sp)); - sp = one.saturating_add(u.div(T::from_num(20)).saturating_mul(sp)); - sp = one.saturating_add(u.div(T::from_num(6)).saturating_mul(sp)); - let sinh_approx = reduced.saturating_mul(sp); - - let mut cp = one; - cp = one.saturating_add(u.div(T::from_num(90)).saturating_mul(cp)); - cp = one.saturating_add(u.div(T::from_num(56)).saturating_mul(cp)); - cp = one.saturating_add(u.div(T::from_num(30)).saturating_mul(cp)); - cp = one.saturating_add(u.div(T::from_num(12)).saturating_mul(cp)); - cp = one.saturating_add(u.div(T::from_num(2)).saturating_mul(cp)); - (sinh_approx, cp) - }; + let mut cp = one; + if T::frac_bits() >= 40 { + sp = one.saturating_add(u.div_int(420).saturating_mul(sp)); + sp = one.saturating_add(u.div_int(342).saturating_mul(sp)); + sp = one.saturating_add(u.div_int(272).saturating_mul(sp)); + sp = one.saturating_add(u.div_int(210).saturating_mul(sp)); + cp = one.saturating_add(u.div_int(462).saturating_mul(cp)); + cp = one.saturating_add(u.div_int(380).saturating_mul(cp)); + cp = one.saturating_add(u.div_int(306).saturating_mul(cp)); + cp = one.saturating_add(u.div_int(240).saturating_mul(cp)); + } + if T::frac_bits() >= 24 { + sp = one.saturating_add(u.div_int(156).saturating_mul(sp)); + sp = one.saturating_add(u.div_int(110).saturating_mul(sp)); + cp = one.saturating_add(u.div_int(182).saturating_mul(cp)); + cp = one.saturating_add(u.div_int(132).saturating_mul(cp)); + } + sp = one.saturating_add(u.div_int(72).saturating_mul(sp)); + sp = one.saturating_add(u.div_int(42).saturating_mul(sp)); + sp = one.saturating_add(u.div_int(20).saturating_mul(sp)); + sp = one.saturating_add(u.div_int(6).saturating_mul(sp)); + cp = one.saturating_add(u.div_int(90).saturating_mul(cp)); + cp = one.saturating_add(u.div_int(56).saturating_mul(cp)); + cp = one.saturating_add(u.div_int(30).saturating_mul(cp)); + cp = one.saturating_add(u.div_int(12).saturating_mul(cp)); + cp = one.saturating_add(u.div_int(2).saturating_mul(cp)); + let (mut sh, mut ch) = (reduced.saturating_mul(sp), cp); // Reconstruct via doubling: sinh(2x) = 2·sinh(x)·cosh(x), // cosh(2x) = cosh²(x) + sinh²(x) @@ -189,18 +184,39 @@ pub fn coth(x: T) -> Result { #[must_use] #[cfg_attr(feature = "verify-no-panic", no_panic::no_panic)] pub fn asinh(x: T) -> T { - if x == T::zero() { - return T::zero(); - } + let zero = T::zero(); + let one = T::one(); - // asinh(x) = atanh(x / sqrt(1 + x²)) - // NonNegative::one_plus_square(x) returns 1 + x², which is always ≥ 1 - let sqrt_term = sqrt_nonneg(NonNegative::one_plus_square(x)); + if x == zero { + return zero; + } - // x / sqrt(1 + x²) is always in (-1, 1) since sqrt(1 + x²) > |x| - let arg = OpenUnitInterval::from_div_by_sqrt_one_plus_square(x, sqrt_term); + let abs_x = x.abs(); + if abs_x <= one { + // asinh(x) = atanh(x / sqrt(1 + x²)); the argument is at most 1/√2, + // inside the CORDIC's direct range, and 1 + x² is always ≥ 1. + let sqrt_term = sqrt_nonneg(NonNegative::one_plus_square(x)); + let arg = OpenUnitInterval::from_div_by_sqrt_one_plus_square(x, sqrt_term); + return atanh_open(arg); + } - atanh_open(arg) + // Beyond 1 that argument approaches 1, where each step of atanh's + // argument reduction triples the error, so use + // asinh(x) = sign(x)·ln(|x|·(1 + sqrt(1 + 1/x²))), + // with 1/x² formed as (1/|x|)², which cannot overflow. + let recip = one.div(abs_x); + let factor = one.saturating_add(sqrt_nonneg(NonNegative::one_plus_square(recip))); + let magnitude = if abs_x > T::max_value() >> 2 { + // |x|·factor would overflow (factor < 2.42): split the logarithm. + ln_positive(abs_x).saturating_add(ln_positive(factor)) + } else { + ln_positive(abs_x.saturating_mul(factor)) + }; + if x.is_negative() { + -magnitude + } else { + magnitude + } } /// Inverse hyperbolic cosine. Domain: `x ≥ 1`. @@ -211,19 +227,34 @@ pub fn asinh(x: T) -> T { #[cfg_attr(feature = "verify-no-panic", no_panic::no_panic)] pub fn acosh(x: T) -> Result { let at_least_one = AtLeastOne::new(x).ok_or_else(|| Error::domain("acosh", "value >= 1"))?; + let zero = T::zero(); + let one = T::one(); - if x == T::one() { - return Ok(T::zero()); + if x == one { + return Ok(zero); } - // acosh(x) = atanh(sqrt(x² - 1) / x) for x > 1 - // NonNegative::square_minus_one gives x² - 1, which is ≥ 0 since x ≥ 1 - let sqrt_term = sqrt_nonneg(NonNegative::square_minus_one(at_least_one)); - - // sqrt(x² - 1) / x is in (-1, 1) for x > 1 since sqrt(x² - 1) < x - let arg = OpenUnitInterval::from_sqrt_square_minus_one_div(sqrt_term, at_least_one); + if x < one.saturating_add(T::half()) { + // acosh(x) = atanh(sqrt(x² - 1) / x); below 1.5 the argument is + // under 0.75, inside the CORDIC's direct range, and x² - 1 ≥ 0. + let sqrt_term = sqrt_nonneg(NonNegative::square_minus_one(at_least_one)); + let arg = OpenUnitInterval::from_sqrt_square_minus_one_div(sqrt_term, at_least_one); + return Ok(atanh_open(arg)); + } - Ok(atanh_open(arg)) + // Beyond 1.5 that argument approaches 1, where each step of atanh's + // argument reduction triples the error, so use + // acosh(x) = ln(x·(1 + sqrt(1 - 1/x²))), + // with 1/x² formed as (1/x)², which cannot overflow. + let recip = UnitInterval::from_reciprocal(at_least_one); + let factor = one.saturating_add(sqrt_nonneg(NonNegative::one_minus_square(recip))); + let result = if x > T::max_value() >> 1 { + // x·factor would overflow (factor < 2): split the logarithm. + ln_positive(x).saturating_add(ln_positive(factor)) + } else { + ln_positive(x.saturating_mul(factor)) + }; + Ok(result) } /// Inverse hyperbolic tangent. Domain: `(-1, 1)`. diff --git a/src/ops/mod.rs b/src/ops/mod.rs index e3cd489..9d7739b 100644 --- a/src/ops/mod.rs +++ b/src/ops/mod.rs @@ -7,7 +7,7 @@ //! //! - [`circular`]: Trigonometric functions (sin, cos, tan, asin, acos, atan, atan2) //! - [`hyperbolic`]: Hyperbolic functions (sinh, cosh, tanh, asinh, acosh, atanh, acoth) -//! - [`exponential`]: Exponential and logarithmic functions (exp, ln, log2, log10, pow2) +//! - [`exponential`]: Exponential and logarithmic functions (exp, ln, log2, log10, pow, pow2) //! - [`algebraic`]: Algebraic functions (sqrt) pub mod algebraic; @@ -18,5 +18,5 @@ pub mod hyperbolic; // Re-export all public functions pub use algebraic::sqrt; pub use circular::{acos, asin, atan, atan2, cos, sin, sin_cos, tan}; -pub use exponential::{exp, ln, log2, log10, pow2}; +pub use exponential::{exp, ln, log2, log10, pow, pow2}; pub use hyperbolic::{acosh, acoth, asinh, atanh, cosh, coth, sinh, sinh_cosh, tanh}; diff --git a/src/tables/chebyshev.rs b/src/tables/chebyshev.rs index 79ca242..8424c30 100644 --- a/src/tables/chebyshev.rs +++ b/src/tables/chebyshev.rs @@ -24,13 +24,18 @@ use crate::traits::CordicNumber; /// /// Stored highest-degree first: `[cₙ, cₙ₋₁, …, c₀]`. /// Evaluates P(x) = cₙ·xⁿ + cₙ₋₁·xⁿ⁻¹ + ··· + c₀ via Horner's method. +/// +/// Uses wrapping arithmetic, so the caller must keep the evaluation in +/// range: for the tables in this module and `0 ≤ x ≤ (π/4)²` every +/// intermediate stays below 1 in magnitude, so wrapping and saturating +/// arithmetic are bit-identical and the overflow checks would be pure cost. #[inline] pub fn horner(coeffs: &[i64; N], x: T) -> T { let mut iter = coeffs.iter(); // First element is the highest-degree coefficient (N ≥ 3 for all tables). let mut result = T::from_i1f63(*iter.next().unwrap_or(&0)); for &coeff in iter { - result = T::from_i1f63(coeff).saturating_add(x.saturating_mul(result)); + result = T::from_i1f63(coeff).wrapping_add(x.wrapping_mul(result)); } result } diff --git a/src/traits.rs b/src/traits.rs index e138679..8a465c3 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -90,23 +90,172 @@ pub trait CordicNumber: /// Saturating subtraction. #[must_use] fn saturating_sub(self, rhs: Self) -> Self; + /// Wrapping multiplication, for use only where overflow is provably + /// impossible: then it matches [`saturating_mul`](Self::saturating_mul) + /// bit for bit without the overflow check. + #[must_use] + fn wrapping_mul(self, rhs: Self) -> Self; + /// Wrapping addition. Same contract as [`wrapping_mul`](Self::wrapping_mul). + #[must_use] + fn wrapping_add(self, rhs: Self) -> Self; + /// Wrapping subtraction. Same contract as [`wrapping_mul`](Self::wrapping_mul). + #[must_use] + fn wrapping_sub(self, rhs: Self) -> Self; + /// Saturating multiplication by an integer. Exact unless it saturates. + #[must_use] + fn mul_int(self, rhs: i32) -> Self; + /// Division by a positive integer, truncated toward zero. Matches + /// [`div`](Self::div) bit for bit for non-negative `self` but divides the + /// raw representation directly, which is several times cheaper on + /// 128-bit types. A zero divisor saturates like [`div`](Self::div). + #[must_use] + fn div_int(self, divisor: u32) -> Self; /// Division. #[must_use] fn div(self, rhs: Self) -> Self; + /// Integer part of the base-2 logarithm, `⌊log₂(self)⌋`, or `None` if + /// `self ≤ 0`. + fn checked_int_log2(self) -> Option; + /// Square root, rounded to the nearest representable value. `self` must + /// be non-negative; the `fixed` types return the root of the magnitude. + #[must_use] + fn sqrt_round(self) -> Self; /// Convert from numeric type. fn from_num(n: N) -> Self; /// Maximum value. fn max_value() -> Self; /// Minimum value. fn min_value() -> Self; - /// Round to nearest integer (half away from zero). + /// Round to nearest integer (half away from zero), saturating at the + /// type's bounds. #[must_use] fn round(self) -> Self; - /// Convert to i32 (truncates toward zero). + /// Convert to i32, rounding toward −∞ and saturating if the value does + /// not fit. #[must_use] fn to_i32(self) -> i32; } +// ============================================================================= +// Square root helpers +// ============================================================================= + +/// Round-to-nearest square root when `X << (FRAC_NBITS + 2)` fits in +/// `$wide`: the raw result is `⌊√N + ½⌋ = (⌊2√N⌋ + 1) >> 1` with +/// `N = X·2^f`, and `⌊2√N⌋ = isqrt(4N)`. +macro_rules! sqrt_round_widened { + ($self:expr, $bits:ty, $wide:ty) => {{ + let n4 = <$wide>::from($self.to_bits().unsigned_abs()) << (Self::FRAC_NBITS + 2); + #[allow( + clippy::cast_possible_truncation, + clippy::cast_possible_wrap, + reason = "the root has at most (total_bits + frac_bits) / 2 + 1 bits" + )] + Self::from_bits(((n4.isqrt() + 1) >> 1) as $bits) + }}; +} + +/// Round-to-nearest square root for 128-bit raw representations, where +/// `X·2^f` does not fit in a machine integer. See [`sqrt_round_u128`]. +macro_rules! sqrt_round_i128 { + ($self:expr, $bits:ty, $wide:ty) => {{ + let root = sqrt_round_u128($self.to_bits().unsigned_abs(), Self::FRAC_NBITS, |a, b| { + // ⌊a·2^f / b⌋ for b ≥ √N is below 2^127, so this cannot + // overflow; the fallback only keeps the function total. + #[allow(clippy::cast_possible_wrap, reason = "operands below 2^127")] + let quotient = + Fixed::checked_div(Self::from_bits(a as i128), Self::from_bits(b as i128)); + #[allow(clippy::cast_sign_loss, reason = "quotient of positive operands")] + { + quotient.map_or(i128::MAX, Self::to_bits) as u128 + } + }); + #[allow(clippy::cast_possible_wrap, reason = "root below 2^127")] + Self::from_bits(root as $bits) + }}; +} + +/// A 256-bit unsigned integer as `(high, low)` limbs. +type U256 = (u128, u128); + +/// `a · b` as a 256-bit product. +const fn wide_mul(a: u128, b: u128) -> U256 { + const MASK: u128 = u64::MAX as u128; + let (a_hi, a_lo) = (a >> 64, a & MASK); + let (b_hi, b_lo) = (b >> 64, b & MASK); + let ll = a_lo * b_lo; + let lh = a_lo * b_hi; + let hl = a_hi * b_lo; + let hh = a_hi * b_hi; + let (mid, mid_carry) = lh.overflowing_add(hl); + let (lo, lo_carry) = ll.overflowing_add(mid << 64); + let hi = hh + (mid >> 64) + ((mid_carry as u128) << 64) + (lo_carry as u128); + (hi, lo) +} + +/// `x << shift` as a 256-bit value, for `0 < shift < 128`. +const fn wide_shl(x: u128, shift: u32) -> U256 { + (x >> (128 - shift), x << shift) +} + +/// `a > b` for 256-bit values. +const fn wide_gt(a: U256, b: U256) -> bool { + a.0 > b.0 || (a.0 == b.0 && a.1 > b.1) +} + +/// Low limb of `a − b`, for `a ≥ b` with a difference below 2^128. +const fn wide_sub_lo(a: U256, b: U256) -> u128 { + a.1.wrapping_sub(b.1) +} + +/// Round-to-nearest square root of `x / 2^frac` (`x < 2^127`), in raw units: +/// `⌊√N + ½⌋` with `N = x·2^frac`, up to 254 bits. +/// +/// `x` is shifted left by the largest `s ≤ leading_zeros(x)` with +/// `s ≡ frac (mod 2)`, giving `q = isqrt(x·2^s)` with 64 significant bits. +/// If `s ≥ frac`, `q >> ((s − frac)/2)` is exactly `⌊√N⌋`. Otherwise +/// `seed = (q + 1) << m` lies above `√N` by less than `2^(m+1)`, one Newton +/// step from above lands on `⌊√N⌋` or `⌊√N⌋ + 1` (its error is below +/// `2^(m − 64)`, `m ≤ 63`), and a 256-bit remainder settles floor and +/// rounding. `div_scaled(a, b)` must return `⌊a·2^frac / b⌋`; it is only +/// called with `b ≥ √N`. +fn sqrt_round_u128(x: u128, frac: u32, div_scaled: impl FnOnce(u128, u128) -> u128) -> u128 { + if x == 0 { + return 0; + } + let lz = x.leading_zeros(); + let shift = lz - ((lz ^ frac) & 1); + let root_hi = (x << shift).isqrt(); + + if shift >= frac { + let excess = (shift - frac) / 2; + if excess >= 1 { + // ⌊2√N⌋ = root_hi >> (excess − 1); round half up. + return ((root_hi >> (excess - 1)) + 1) >> 1; + } + // root_hi = ⌊√N⌋, and N fits in 128 bits since shift = frac ≤ lz. + let root_sq = root_hi * root_hi; + let rem = (x << frac) - root_sq; + return if rem > root_hi { root_hi + 1 } else { root_hi }; + } + + let missing = (frac - shift) / 2; + #[allow(clippy::cast_sign_loss, reason = "i128::MAX is positive")] + let seed = ((root_hi + 1) << missing).min(i128::MAX as u128); + let newton = u128::midpoint(seed, div_scaled(x, seed)); + + let n = wide_shl(x, frac); + let newton_sq = wide_mul(newton, newton); + let (root, rem) = if wide_gt(newton_sq, n) { + // newton = ⌊√N⌋ + 1: N − (newton − 1)² = (2·newton − 1) − (newton² − N). + (newton - 1, (2 * newton - 1) - wide_sub_lo(newton_sq, n)) + } else { + (newton, wide_sub_lo(n, newton_sq)) + }; + // rem = N − root² ∈ [0, 2·root]; round up exactly when rem > root. + if rem > root { root + 1 } else { root } +} + // ============================================================================= // Generic implementations using macros // ============================================================================= @@ -122,6 +271,8 @@ macro_rules! impl_cordic_generic { ( $fixed_type:ident, $bits_type:ty, + $wide_type:ty, // Unsigned type holding bits << (frac + 2) for sqrt + $sqrt_impl:ident, // sqrt_round_widened or sqrt_round_i128 $total_bits:expr, $max_frac:ty, // Maximum fractional bits for the type $pi_frac:ty, // Max frac bits where PI fits (total - 2) @@ -240,6 +391,59 @@ macro_rules! impl_cordic_generic { Fixed::saturating_sub(self, rhs) } + #[inline] + fn wrapping_mul(self, rhs: Self) -> Self { + Fixed::wrapping_mul(self, rhs) + } + + #[inline] + fn wrapping_add(self, rhs: Self) -> Self { + Fixed::wrapping_add(self, rhs) + } + + #[inline] + fn wrapping_sub(self, rhs: Self) -> Self { + Fixed::wrapping_sub(self, rhs) + } + + #[inline] + fn mul_int(self, rhs: i32) -> Self { + match <$bits_type>::try_from(rhs) { + Ok(k) => Fixed::saturating_mul_int(self, k), + // |rhs| exceeds the raw type's range (8- and 16-bit types + // only), so the exact product saturates unless self is 0. + Err(_) => { + if self == Self::ZERO { + Self::ZERO + } else if self.is_negative() != (rhs < 0) { + Self::MIN + } else { + Self::MAX + } + } + } + } + + #[inline] + fn div_int(self, divisor: u32) -> Self { + match <$bits_type>::try_from(divisor) { + Ok(k) => match Fixed::checked_div_int(self, k) { + Some(v) => v, + // Division by zero: saturate based on sign. + None => { + if self.is_negative() { + Self::MIN + } else { + Self::MAX + } + } + }, + // The divisor exceeds the raw type's range, so it exceeds + // |self| in raw units and the truncated quotient is zero. + Err(_) => Self::ZERO, + } + } + #[inline] fn div(self, rhs: Self) -> Self { match Fixed::checked_div(self, rhs) { @@ -255,6 +459,16 @@ macro_rules! impl_cordic_generic { } } + #[inline] + fn checked_int_log2(self) -> Option { + Fixed::checked_int_log2(self) + } + + #[inline] + fn sqrt_round(self) -> Self { + $sqrt_impl!(self, $bits_type, $wide_type) + } + #[inline] fn from_num(n: N) -> Self { Self::from_num(n) @@ -272,16 +486,12 @@ macro_rules! impl_cordic_generic { #[inline] fn round(self) -> Self { - Fixed::round(self) + Fixed::saturating_round(self) } #[inline] - #[allow( - clippy::cast_possible_truncation, - reason = "intentional truncation to target type" - )] fn to_i32(self) -> i32 { - self.to_num::() + self.saturating_to_num::() } } }; @@ -298,28 +508,69 @@ use fixed::types::extra::{ // - For PI (~3.14), need 2 integer bits, so Fract ≤ 6 (I2F6) // - For FRAC_PI_2, FRAC_PI_4, LN_2, need 1 integer bit, so Fract ≤ 7 (I1F7) // Being conservative: require Fract ≤ 5 so we have headroom -impl_cordic_generic!(FixedI8, i8, 8, U8, U5, U6, U7); +impl_cordic_generic!(FixedI8, i8, u16, sqrt_round_widened, 8, U8, U5, U6, U7); // FixedI16: 16 total bits // - For PI, need Fract ≤ 14 (I2F14) // - For FRAC_PI_2, FRAC_PI_4, LN_2, need Fract ≤ 15 (I1F15) // - Conservative: Fract ≤ 13 -impl_cordic_generic!(FixedI16, i16, 16, U16, U13, U14, U15); +impl_cordic_generic!( + FixedI16, + i16, + u32, + sqrt_round_widened, + 16, + U16, + U13, + U14, + U15 +); // FixedI32: 32 total bits // - For PI, need Fract ≤ 30 // - For FRAC_PI_2, FRAC_PI_4, LN_2, need Fract ≤ 31 // - Conservative: Fract ≤ 29 -impl_cordic_generic!(FixedI32, i32, 32, U32, U29, U30, U31); +impl_cordic_generic!( + FixedI32, + i32, + u64, + sqrt_round_widened, + 32, + U32, + U29, + U30, + U31 +); // FixedI64: 64 total bits // - For PI, need Fract ≤ 62 // - For FRAC_PI_2, FRAC_PI_4, LN_2, need Fract ≤ 63 // - Conservative: Fract ≤ 61 -impl_cordic_generic!(FixedI64, i64, 64, U64, U61, U62, U63); +impl_cordic_generic!( + FixedI64, + i64, + u128, + sqrt_round_widened, + 64, + U64, + U61, + U62, + U63 +); // FixedI128: 128 total bits // - For PI, need Fract ≤ 126 // - For FRAC_PI_2, FRAC_PI_4, LN_2, need Fract ≤ 127 // - Conservative: Fract ≤ 125 -impl_cordic_generic!(FixedI128, i128, 128, U128, U125, U126, U127); +// No wider machine integer exists, so sqrt takes the Newton path. +impl_cordic_generic!( + FixedI128, + i128, + u128, + sqrt_round_i128, + 128, + U128, + U125, + U126, + U127 +); diff --git a/tests/unit/kernel/cordic.rs b/tests/unit/kernel/cordic.rs index 617901d..bb0d56d 100644 --- a/tests/unit/kernel/cordic.rs +++ b/tests/unit/kernel/cordic.rs @@ -27,3 +27,54 @@ mod tests { ); } } + +/// Sweep both early-terminating kernels against `f64` at three widths. +#[cfg(test)] +mod early_termination { + use crate::unit::support::Lcg; + use fixed::types::{I16F16, I32F32, I64F64}; + use fixed_analytics::CordicNumber; + use fixed_analytics::kernel::{circular_vectoring, hyperbolic_vectoring}; + + fn sweep(tol_circ: f64, tol_hyp: f64) { + let mut rng = Lcg(0xC0DE); + for i in 0..=2000 { + let v = if i <= 1000 { + -1.0 + 2.0 * f64::from(i) / 1000.0 + } else { + rng.range(-1.0, 1.0) + }; + let y = ::from_num(v); + let z: f64 = circular_vectoring(T::one(), y, T::zero()).2.to_num(); + let want = v.atan(); + assert!( + (z - want).abs() < tol_circ, + "atan({v}) via kernel = {z}, want {want} (tol {tol_circ})" + ); + if v.abs() <= 0.75 { + let zh: f64 = hyperbolic_vectoring(T::one(), y, T::zero()).2.to_num(); + let want_h = v.atanh(); + assert!( + (zh - want_h).abs() < tol_hyp, + "atanh({v}) via kernel = {zh}, want {want_h} (tol {tol_hyp})" + ); + } + } + } + + #[test] + fn kernels_track_f64_at_i16f16() { + sweep::(8e-5, 8e-5); + } + + #[test] + fn kernels_track_f64_at_i32f32() { + sweep::(3e-9, 3e-9); + } + + #[test] + fn kernels_track_f64_at_i64f64() { + // Above an f64 ulp: the libm reference varies by platform. + sweep::(6e-16, 6e-16); + } +} diff --git a/tests/unit/mod.rs b/tests/unit/mod.rs index 0ee5e5f..f360249 100644 --- a/tests/unit/mod.rs +++ b/tests/unit/mod.rs @@ -4,6 +4,7 @@ mod error; mod kernel; mod ops; mod smoke; +mod support; mod tables; mod traits; mod verification; diff --git a/tests/unit/ops/algebraic.rs b/tests/unit/ops/algebraic.rs index 1813d0e..fa032b5 100644 --- a/tests/unit/ops/algebraic.rs +++ b/tests/unit/ops/algebraic.rs @@ -56,3 +56,217 @@ mod tests { } } } + +/// Exactness: every result is the representable value nearest the true +/// root, checked against the `fixed` crate's rounded-down square root plus +/// a 256-bit remainder test for the rounding direction. +#[cfg(test)] +#[allow( + clippy::unwrap_used, + clippy::cast_sign_loss, + clippy::cast_possible_wrap, + clippy::cast_possible_truncation, + reason = "test code" +)] +mod exact { + use crate::unit::support::Lcg; + use fixed::traits::Fixed; + use fixed::types::{ + I3F125, I4F4, I4F60, I8F8, I8F24, I16F16, I24F8, I32F32, I48F16, I64F64, I128F0, + }; + use fixed_analytics::sqrt; + + type U256 = (u128, u128); + + fn wide_mul(a: u128, b: u128) -> U256 { + const MASK: u128 = u64::MAX as u128; + let (a_hi, a_lo) = (a >> 64, a & MASK); + let (b_hi, b_lo) = (b >> 64, b & MASK); + let (mid, mid_carry) = (a_lo * b_hi).overflowing_add(a_hi * b_lo); + let (lo, lo_carry) = (a_lo * b_lo).overflowing_add(mid << 64); + let hi = a_hi * b_hi + (mid >> 64) + (u128::from(mid_carry) << 64) + u128::from(lo_carry); + (hi, lo) + } + + fn wide_shl(x: u128, shift: u32) -> U256 { + if shift == 0 { + (0, x) + } else { + (x >> (128 - shift), x << shift) + } + } + + /// `a - b` for `a ≥ b`. + fn wide_sub(a: U256, b: U256) -> U256 { + let (lo, borrow) = a.1.overflowing_sub(b.1); + (a.0 - b.0 - u128::from(borrow), lo) + } + + /// Raw bits of the correctly rounded root of `x`. + fn expected_bits(x: T) -> i128 + where + T::Bits: Into, + { + let floor: i128 = x.sqrt().to_bits().into(); + let n = wide_shl(x.to_bits().into() as u128, T::FRAC_NBITS); + let f = floor as u128; + let rem = wide_sub(n, wide_mul(f, f)); + assert_eq!(rem.0, 0, "remainder exceeds 2·root + 1"); + if rem.1 > f { floor + 1 } else { floor } + } + + /// Deterministic sample set spanning the whole non-negative range of `T`: + /// special values, powers of two and neighbours, near-squares, and + /// uniformly random raw bits at several magnitudes. + fn samples(random: usize) -> Vec + where + T::Bits: TryFrom + Into, + { + let bits = T::INT_NBITS + T::FRAC_NBITS; + let max: i128 = T::MAX.to_bits().into(); + let mut raw: Vec = vec![0, 1, 2, 3, max, max - 1]; + for e in 0..(bits - 1) { + for d in [-2, -1, 0, 1, 2] { + raw.push((1i128 << e) + d); + } + } + // X = k² ± j lands √(X·2^f) near an integer, exercising both + // rounding directions and the 128-bit Newton correction step. + for p in 1..64u32 { + let k = 1i128 << p; + if let Some(k2) = k.checked_mul(k) { + for j in [-3, -2, -1, 0, 1, 2, 3] { + raw.push(k2 + j); + } + } + } + let mut rng = Lcg(0xABCD_EF01 ^ u64::from(bits)); + for _ in 0..random { + let r = (u128::from(rng.next_u64()) | (u128::from(rng.next_u64()) << 64)) as i128; + let full = (r >> (128 - bits)).abs().min(max); + raw.extend([ + full, + full >> (bits / 2), + full >> (bits / 3), + full >> (2 * bits / 3), + ]); + let k = full >> (bits / 2 + 1); + for j in [-1, 0, 1] { + raw.push(k * k + k + j); + } + } + raw.into_iter() + .filter(|&v| (0..=max).contains(&v)) + .filter_map(|v| T::Bits::try_from(v).ok()) + .map(T::from_bits) + .collect() + } + + fn check(random: usize) + where + T: Fixed + fixed_analytics::CordicNumber + core::fmt::Display, + T::Bits: TryFrom + Into, + { + for x in samples::(random) { + let got: i128 = sqrt(x).unwrap().to_bits().into(); + let want = expected_bits(x); + assert_eq!(got, want, "sqrt({x}) raw bits: got {got}, want {want}"); + } + } + + #[test] + fn i4f4_and_i8f8_are_correctly_rounded() { + check::(500); + check::(2000); + } + + #[test] + fn i16f16_i8f24_i24f8_are_correctly_rounded() { + check::(5000); + check::(2000); + check::(2000); + } + + #[test] + fn i32f32_i4f60_i48f16_are_correctly_rounded() { + check::(5000); + check::(2000); + check::(2000); + } + + #[test] + fn i64f64_i3f125_i128f0_are_correctly_rounded() { + // 128-bit types take the Newton path above 2^(128 - 2·frac). + check::(5000); + check::(3000); + check::(1000); + } + + /// Below 0.25 on I64F64 the exact root fits in `u128` arithmetic. + #[test] + fn i64f64_sub_unit_matches_integer_root() { + let mut rng = Lcg(0x5EED); + let mut raw: Vec = vec![1, 2, 3, u64::MAX.into()]; + // Log-spaced from one ulp (2^-64 ≈ 5e-20) up to 0.25. + for i in 0..4000 { + let exponent = 62.0 * f64::from(i) / 4000.0; + raw.push((exponent).exp2().max(1.0) as u128); + raw.push(rng.next_u64().into()); + raw.push(u128::from(rng.next_u64()) >> (rng.next_u64() % 62)); + } + for x_raw in raw { + let x_raw = x_raw.min((1u128 << 62) - 1); + let x = I64F64::from_bits(x_raw as i128); + // round(√(X·2^64)) = (⌊2·√(X·2^64)⌋ + 1) >> 1 = (isqrt(X·2^66) + 1) >> 1 + let want = (((x_raw << 66).isqrt() + 1) >> 1) as i128; + let got = sqrt(x).unwrap().to_bits(); + assert_eq!(got, want, "sqrt({x}) raw bits: got {got}, want {want}"); + } + } + + /// Inputs just below a perfect square make the Newton step land one + /// above the floor. + #[test] + fn i64f64_newton_correction_cases() { + for p in 33..63u32 { + for j in 1..=3i128 { + let x = I64F64::from_bits((1i128 << (2 * p - 64)) - j); + let got: i128 = sqrt(x).unwrap().to_bits(); + let want = expected_bits(x); + assert_eq!(got, want, "sqrt({x}) raw bits: got {got}, want {want}"); + } + } + } + + #[test] + fn special_values() { + assert_eq!(sqrt(I64F64::ZERO).unwrap(), I64F64::ZERO); + assert_eq!(sqrt(I64F64::ONE).unwrap(), I64F64::ONE); + assert_eq!(sqrt(I64F64::from_num(4)).unwrap(), I64F64::from_num(2)); + assert_eq!(sqrt(I64F64::from_num(0.25)).unwrap(), I64F64::from_num(0.5)); + assert_eq!( + sqrt(I64F64::from_bits(1)).unwrap(), + I64F64::from_bits(1 << 32) + ); + assert_eq!( + sqrt(I64F64::MAX).unwrap().to_bits(), + expected_bits(I64F64::MAX) + ); + assert_eq!( + sqrt(I3F125::MAX).unwrap().to_bits(), + expected_bits(I3F125::MAX) + ); + assert_eq!( + i128::from(sqrt(I16F16::MAX).unwrap().to_bits()), + expected_bits(I16F16::MAX) + ); + } + + #[test] + fn asin_acos_still_saturate_near_one() { + let x = I64F64::ONE - I64F64::from_bits(1); + assert!( + (fixed_analytics::asin(x).unwrap() - I64F64::FRAC_PI_2).abs() < I64F64::from_num(1e-9) + ); + } +} diff --git a/tests/unit/ops/circular.rs b/tests/unit/ops/circular.rs index 9de8625..4338348 100644 --- a/tests/unit/ops/circular.rs +++ b/tests/unit/ops/circular.rs @@ -516,3 +516,102 @@ mod tests { } } } + +/// Angle reduction near the type bounds and inverse-function accuracy. +#[cfg(test)] +#[allow(clippy::unwrap_used, reason = "test code uses unwrap for conciseness")] +mod reduction_and_accuracy { + use crate::unit::support::Lcg; + use fixed::types::{I24F8, I32F32, I64F64}; + use fixed_analytics::{CordicNumber, acos, asin, atan, atan2, sin_cos}; + + /// Reference reduction with the type's own 2π, in f64. + fn reduced_f64(angle: T) -> f64 { + let two_pi: f64 = (T::pi() + T::pi()).to_num(); + let a: f64 = angle.to_num(); + (a / two_pi).round().mul_add(-two_pi, a) + } + + #[test] + fn sin_cos_is_correct_within_pi_of_max() { + // n·2π exceeds MAX here; saturation used to return (0, 1). + for angle in [I32F32::MAX, I32F32::MIN, I32F32::MAX - I32F32::ONE] { + let (s, c) = sin_cos(angle); + let r = reduced_f64(angle); + let (s, c): (f64, f64) = (s.to_num(), c.to_num()); + assert!( + (s - r.sin()).abs() < 1e-5, + "sin({angle}) = {s}, want {}", + r.sin() + ); + assert!( + (c - r.cos()).abs() < 1e-5, + "cos({angle}) = {c}, want {}", + r.cos() + ); + } + for angle in [I24F8::MAX, I24F8::MIN] { + let (s, c) = sin_cos(angle); + let r = reduced_f64(angle); + let (s, c): (f64, f64) = (s.to_num(), c.to_num()); + assert!( + (s - r.sin()).abs() < 0.02, + "sin({angle}) = {s}, want {}", + r.sin() + ); + assert!( + (c - r.cos()).abs() < 0.02, + "cos({angle}) = {c}, want {}", + r.cos() + ); + } + } + + #[test] + fn sin_cos_large_angles_i64f64() { + let mut rng = Lcg(0x51); + for _ in 0..500 { + let v = rng.range(-1e6, 1e6); + let angle = I64F64::from_num(v); + let (s, c) = sin_cos(angle); + let r = reduced_f64(angle); + let (s, c): (f64, f64) = (s.to_num(), c.to_num()); + assert!( + (s - r.sin()).abs() < 1e-9, + "sin({v}) = {s}, want {}", + r.sin() + ); + assert!( + (c - r.cos()).abs() < 1e-9, + "cos({v}) = {c}, want {}", + r.cos() + ); + } + } + + #[test] + fn inverse_functions_track_f64_at_i64f64() { + let mut rng = Lcg(0xA7A); + for _ in 0..2000 { + let v = rng.range(-1.0, 1.0); + let x = I64F64::from_num(v); + let asin_v: f64 = asin(x).unwrap().to_num(); + assert!((asin_v - v.asin()).abs() < 2e-15, "asin({v}) = {asin_v}"); + let acos_v: f64 = acos(x).unwrap().to_num(); + assert!((acos_v - v.acos()).abs() < 2e-15, "acos({v}) = {acos_v}"); + let big = (rng.range(-30.0, 30.0)).exp2() * v.signum(); + let atan_v: f64 = atan(I64F64::from_num(big)).to_num(); + assert!( + (atan_v - big.atan()).abs() < 1e-15, + "atan({big}) = {atan_v}" + ); + let y = rng.range(-5.0, 5.0); + let x2 = rng.range(-5.0, 5.0); + let atan2_v: f64 = atan2(I64F64::from_num(y), I64F64::from_num(x2)).to_num(); + assert!( + (atan2_v - y.atan2(x2)).abs() < 1e-15, + "atan2({y}, {x2}) = {atan2_v}" + ); + } + } +} diff --git a/tests/unit/ops/exponential.rs b/tests/unit/ops/exponential.rs index b62a74a..d220795 100644 --- a/tests/unit/ops/exponential.rs +++ b/tests/unit/ops/exponential.rs @@ -424,3 +424,199 @@ mod tests { } } } + +/// Argument reduction, saturation, and the bit-scan normalisation in `ln`. +#[cfg(test)] +#[allow(clippy::unwrap_used, reason = "test code uses unwrap for conciseness")] +mod reduction { + use crate::unit::support::Lcg; + use fixed::types::{I4F60, I16F16, I24F8, I32F32, I64F64}; + use fixed_analytics::{CordicNumber, exp, ln, log2, pow, pow2, sqrt}; + + #[test] + fn exp_saturates_at_the_extremes_of_wide_types() { + // x / ln2 does not fit in i32 here; the conversion saturates. + assert_eq!(exp(I32F32::MAX), I32F32::MAX); + assert_eq!(exp(I32F32::MIN), I32F32::ZERO); + assert_eq!(exp(I32F32::from_num(2e9)), I32F32::MAX); + assert_eq!(exp(I64F64::MAX), I64F64::MAX); + assert_eq!(exp(I64F64::MIN), I64F64::ZERO); + assert_eq!(exp(I64F64::from_num(1.5e9)), I64F64::MAX); + assert_eq!(exp(I64F64::from_num(-1.5e9)), I64F64::ZERO); + assert_eq!(exp(I16F16::MAX), I16F16::MAX); + assert_eq!(exp(I16F16::MIN), I16F16::ZERO); + assert_eq!(pow2(I64F64::MAX), I64F64::MAX); + assert_eq!(pow2(I64F64::MIN), I64F64::ZERO); + } + + #[test] + fn exp_thresholds_after_correction_i16f16() { + // k reaches 15 (overflow) or -18 (zero) only after the correction. + assert_eq!(exp(I16F16::from_num(10.5)), I16F16::MAX); + assert!(exp(I16F16::from_num(10.39)) < I16F16::MAX); + assert_eq!(exp(I16F16::from_num(-12.1)), I16F16::ZERO); + assert!(exp(I16F16::from_num(-11.4)) > I16F16::ZERO); + } + + #[test] + fn exp_estimate_corrected_downward() { + // Just below a negative multiple of ln2 the reciprocal estimate + // rounds up to the multiple; r < 0 then steps k down. + let x16 = I16F16::LN_2.mul_int(-3) - I16F16::from_bits(1); + let got16: f64 = exp(x16).to_num(); + assert!((got16 - 0.125).abs() < 3e-5, "exp({x16}) = {got16}"); + let x64 = I64F64::LN_2.mul_int(-3) - I64F64::from_bits(1); + let got64: f64 = exp(x64).to_num(); + assert!((got64 - 0.125).abs() < 1e-12, "exp({x64}) = {got64}"); + } + + #[test] + fn exp_tracks_f64_at_i64f64() { + let mut rng = Lcg(0xE4); + for i in 0..2000 { + let v = if i < 1000 { + -40.0 + 80.0 * f64::from(i) / 1000.0 + } else { + rng.range(-8.0, 8.0) + }; + let x = I64F64::from_num(v); + let got: f64 = exp(x).to_num(); + let want = x.to_num::().exp(); + // The f64 reference loses 2^-53·|x| of relative precision in x. + let tol = 8.0f64.mul_add(2f64.powi(-64), 5e-16 * (1.0 + v.abs()) * want); + assert!((got - want).abs() < tol, "exp({v}) = {got}, want {want}"); + } + } + + #[test] + fn exp_works_on_types_with_few_integer_bits() { + // I4F60 (range ±8) cannot hold the Taylor divisors or k·ln2 in T. + for i in 0..=200 { + let v = -2.0 + 4.0 * f64::from(i) / 200.0; + let x = I4F60::from_num(v); + let got: f64 = exp(x).to_num(); + let want = x.to_num::().exp(); + assert!( + ((got - want) / want).abs() < 5e-15, + "I4F60 exp({v}) = {got}, want {want}" + ); + } + assert_eq!(exp(I4F60::from_num(3)), I4F60::MAX); + } + + #[test] + fn exp_on_a_type_with_more_integer_than_fractional_bits() { + // I24F8 divides for the estimate; k·ln2 carries k·0.0017 of quantisation. + for i in 0..=100 { + let v = -5.0 + 8.0 * f64::from(i) / 100.0; + let x = I24F8::from_num(v); + let got: f64 = exp(x).to_num(); + let want = x.to_num::().exp(); + assert!( + (got - want).abs() < 0.01 * want.max(1.0), + "I24F8 exp({v}) = {got}, want {want}" + ); + } + assert_eq!(exp(I24F8::MAX), I24F8::MAX); + assert_eq!(exp(I24F8::MIN), I24F8::ZERO); + } + + #[test] + fn ln_normalisation_branches() { + let ulp = I64F64::from_bits(1); + let cases: [(I64F64, &str); 8] = [ + (I64F64::from_num(2), "exactly two"), + (I64F64::from_num(4) + ulp, "lands on two after one shift"), + ( + I64F64::from_num(4) + ulp + ulp, + "lands on 1 + ulp after two shifts", + ), + (I64F64::from_num(5), "in (2, 4)"), + (I64F64::from_num(1000), "large"), + (I64F64::from_num(0.3), "below one half"), + (I64F64::from_num(0.25), "exactly one quarter"), + (I64F64::from_bits(1), "one ulp"), + ]; + for (x, label) in cases { + let got: f64 = ln(x).unwrap().to_num(); + let want: f64 = x.to_num::().ln(); + assert!( + (got - want).abs() < 1e-15 * want.abs().max(1.0), + "ln({x}) [{label}] = {got}, want {want}" + ); + } + for k in 1..=60 { + let want = f64::from(k) * core::f64::consts::LN_2; + let got_log2: f64 = log2(I64F64::from_num(1u64 << k)).unwrap().to_num(); + assert!( + (got_log2 - f64::from(k)).abs() < 1e-14, + "log2(2^{k}) = {got_log2}" + ); + let got_ln: f64 = ln(I64F64::from_num(1u64 << k)).unwrap().to_num(); + assert!( + (got_ln - want).abs() < 1e-14, + "ln(2^{k}) = {got_ln}, want {want}" + ); + } + assert!(ln(I16F16::from_num(4) + I16F16::from_bits(1)).is_ok()); + } + + #[test] + fn pow_special_cases_and_domain() { + assert_eq!(pow(I16F16::ZERO, I16F16::ZERO).unwrap(), I16F16::ONE); + assert_eq!(pow(I16F16::from_num(7), I16F16::ZERO).unwrap(), I16F16::ONE); + assert_eq!( + pow(I16F16::ZERO, I16F16::from_num(2)).unwrap(), + I16F16::ZERO + ); + assert_eq!(pow(I16F16::ONE, I16F16::from_num(-9)).unwrap(), I16F16::ONE); + assert!(pow(I16F16::ZERO, -I16F16::ONE).is_err()); + assert!(pow(-I16F16::ONE, I16F16::from_num(2)).is_err()); + assert_eq!( + pow(I16F16::from_num(200), I16F16::from_num(3)).unwrap(), + I16F16::MAX + ); + let got: f64 = pow(I16F16::from_num(2), I16F16::from_num(10)) + .unwrap() + .to_num(); + assert!((got - 1024.0).abs() < 0.5, "I16F16 2^10 = {got}"); + } + + #[test] + fn pow_tracks_f64_at_i64f64() { + let mut rng = Lcg(0x90); + for _ in 0..2000 { + let base = I64F64::from_num((rng.range(-20.0, 20.0)).exp2()); + let exponent = I64F64::from_num(rng.range(-3.0, 3.0)); + let (b, e): (f64, f64) = (base.to_num(), exponent.to_num()); + let got: f64 = pow(base, exponent).unwrap().to_num(); + let want = b.powf(e); + // Small results are quantised to a few raw ulps. + let tol = 8.0f64.mul_add(2f64.powi(-64), 1e-14 * want); + assert!( + (got - want).abs() < tol, + "pow({b}, {e}) = {got}, want {want}" + ); + } + let x = I64F64::from_num(2.5); + let root: f64 = pow(x, I64F64::from_num(0.5)).unwrap().to_num(); + let want: f64 = sqrt(x).unwrap().to_num(); + assert!((root - want).abs() < 1e-15, "2.5^0.5 = {root}, want {want}"); + let inv: f64 = pow(x, -I64F64::ONE).unwrap().to_num(); + assert!((inv - 0.4).abs() < 1e-15, "2.5^-1 = {inv}"); + } + + #[test] + fn ln_tracks_f64_at_i64f64() { + let mut rng = Lcg(0x1A); + for _ in 0..3000 { + let x = I64F64::from_num((rng.range(-60.0, 62.0)).exp2()); + let got: f64 = ln(x).unwrap().to_num(); + let want = x.to_num::().ln(); + assert!( + (got - want).abs() < 1e-15 * want.abs().max(1.0), + "ln({x}) = {got}, want {want}" + ); + } + } +} diff --git a/tests/unit/ops/hyperbolic.rs b/tests/unit/ops/hyperbolic.rs index 977f1ad..0dc0f65 100644 --- a/tests/unit/ops/hyperbolic.rs +++ b/tests/unit/ops/hyperbolic.rs @@ -487,3 +487,113 @@ mod tests { } } } + +/// Integer-bit-poor types, and the logarithmic inverse forms. +#[cfg(test)] +#[allow(clippy::unwrap_used, reason = "test code uses unwrap for conciseness")] +mod wide_and_narrow { + use crate::unit::support::Lcg; + use fixed::types::{I4F12, I4F60, I16F16, I32F32, I64F64}; + use fixed_analytics::{acosh, asinh, sinh_cosh}; + + #[test] + fn sinh_cosh_works_on_types_with_few_integer_bits() { + // I4F12 and I4F60 (range ±8) cannot hold the Taylor divisors in T. + for i in 0..=100 { + let v = -1.5 + 3.0 * f64::from(i) / 100.0; + let x12 = I4F12::from_num(v); + let (s12, c12) = sinh_cosh(x12); + let (s12, c12, v12): (f64, f64, f64) = (s12.to_num(), c12.to_num(), x12.to_num()); + assert!((s12 - v12.sinh()).abs() < 5e-3, "I4F12 sinh({v12}) = {s12}"); + assert!((c12 - v12.cosh()).abs() < 5e-3, "I4F12 cosh({v12}) = {c12}"); + let x60 = I4F60::from_num(v); + let (s60, c60) = sinh_cosh(x60); + let (s60, c60, v60): (f64, f64, f64) = (s60.to_num(), c60.to_num(), x60.to_num()); + assert!( + (s60 - v60.sinh()).abs() < 4e-15, + "I4F60 sinh({v60}) = {s60}" + ); + assert!( + (c60 - v60.cosh()).abs() < 4e-15, + "I4F60 cosh({v60}) = {c60}" + ); + } + } + + #[test] + fn asinh_acosh_track_f64_at_i64f64() { + let mut rng = Lcg(0xA5); + for _ in 0..2000 { + let sign = if rng.unit() < 0.5 { -1.0 } else { 1.0 }; + let x = I64F64::from_num(sign * (rng.range(-20.0, 60.0)).exp2()); + let got: f64 = asinh(x).to_num(); + let want = x.to_num::().asinh(); + assert!( + (got - want).abs() < 1e-15 * want.abs().max(1.0), + "asinh({x}) = {got}, want {want}" + ); + // acosh is ill-conditioned near 1; start where f64 is trustworthy. + let w = I64F64::from_num(1.0 + (rng.range(-10.0, 61.0)).exp2()); + let got_w: f64 = acosh(w).unwrap().to_num(); + let want_w = w.to_num::().acosh(); + assert!( + (got_w - want_w).abs() < 1e-13 * want_w.abs().max(1.0), + "acosh({w}) = {got_w}, want {want_w}" + ); + } + // Both branches around the 1.5 threshold of acosh. + for w in [1.01f64, 1.25, 1.49, 1.5, 1.51, 2.0] { + let got: f64 = acosh(I64F64::from_num(w)).unwrap().to_num(); + assert!((got - w.acosh()).abs() < 1e-9, "acosh({w}) = {got}"); + } + } + + #[test] + fn asinh_acosh_near_the_type_bounds() { + // The product would overflow, so the logarithm is split. + for (x, want) in [ + (I16F16::MAX, (2.0 * I16F16::MAX.to_num::()).ln()), + (I16F16::from_num(10_000), 10_000f64.asinh()), + (I16F16::MIN, -(2.0 * I16F16::MAX.to_num::()).ln()), + ] { + let got: f64 = asinh(x).to_num(); + assert!((got - want).abs() < 2e-4, "asinh({x}) = {got}, want {want}"); + } + for x in [I16F16::MAX, I16F16::from_num(20_000)] { + let want = (2.0f64 * x.to_num::()).ln(); + let got: f64 = acosh(x).unwrap().to_num(); + assert!((got - want).abs() < 2e-4, "acosh({x}) = {got}, want {want}"); + } + let want_32 = (2.0f64 * I32F32::MAX.to_num::()).ln(); + let got_32: f64 = acosh(I32F32::MAX).unwrap().to_num(); + assert!( + (got_32 - want_32).abs() < 1e-8, + "acosh(I32F32::MAX) = {got_32}" + ); + let want_64 = (2.0f64 * I64F64::MAX.to_num::()).ln(); + let got_64: f64 = acosh(I64F64::MAX).unwrap().to_num(); + assert!( + (got_64 - want_64).abs() < 1e-12, + "acosh(I64F64::MAX) = {got_64}" + ); + let got_asinh: f64 = asinh(I64F64::MAX).to_num(); + assert!( + (got_asinh - want_64).abs() < 1e-12, + "asinh(MAX) = {got_asinh}, want {want_64}" + ); + assert_eq!(asinh(-I64F64::from_num(3)), -asinh(I64F64::from_num(3))); + } + + #[test] + fn asinh_is_accurate_for_moderate_arguments_at_i16f16() { + for i in 1..=200 { + let v = f64::from(i).mul_add(0.1, 1.0); + let got: f64 = asinh(I16F16::from_num(v)).to_num(); + let want = v.asinh(); + assert!( + (got - want).abs() < 16.0 / 65_536.0, + "asinh({v}) = {got}, want {want}" + ); + } + } +} diff --git a/tests/unit/support.rs b/tests/unit/support.rs new file mode 100644 index 0000000..9ad7257 --- /dev/null +++ b/tests/unit/support.rs @@ -0,0 +1,29 @@ +//! Shared helpers for the unit tests. + +/// Deterministic generator for reproducible sweeps. +pub struct Lcg(pub u64); + +impl Lcg { + pub const fn next_u64(&mut self) -> u64 { + self.0 = self + .0 + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + let x = self.0; + (x ^ (x >> 29)).wrapping_mul(0x9E37_79B9_7F4A_7C15) ^ (x >> 32) + } + + /// Uniform in [0, 1). + #[allow( + clippy::cast_precision_loss, + reason = "53 bits of randomness is plenty" + )] + pub const fn unit(&mut self) -> f64 { + (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64 + } + + /// Uniform in [lo, hi). + pub fn range(&mut self, lo: f64, hi: f64) -> f64 { + (hi - lo).mul_add(self.unit(), lo) + } +} diff --git a/tests/unit/tables/chebyshev.rs b/tests/unit/tables/chebyshev.rs index efe2d9a..566826d 100644 --- a/tests/unit/tables/chebyshev.rs +++ b/tests/unit/tables/chebyshev.rs @@ -73,3 +73,70 @@ mod tests { assert_eq!(COS_Q_HI.len(), 7); } } + +/// `horner` wraps on the strength of every intermediate staying below 1 in +/// magnitude: check the bound and the bit-identity with saturating Horner. +#[cfg(test)] +#[allow( + clippy::cast_precision_loss, + clippy::unwrap_used, + reason = "test code uses f64 for verification" +)] +mod wrapping_invariant { + use fixed::types::{I16F16, I32F32, I64F64}; + use fixed_analytics::CordicNumber; + use fixed_analytics::tables::chebyshev::{COS_Q_HI, COS_Q_LO, SIN_P_HI, SIN_P_LO, horner}; + + const SCALE: f64 = (1_u64 << 63) as f64; + const U_MAX: f64 = core::f64::consts::FRAC_PI_4 * core::f64::consts::FRAC_PI_4; + + fn max_intermediate(coeffs: &[i64]) -> f64 { + let mut worst: f64 = 0.0; + for i in 0..=1000 { + let u = U_MAX * f64::from(i) / 1000.0; + let mut iter = coeffs.iter(); + let mut acc = *iter.next().unwrap_or(&0) as f64 / SCALE; + worst = worst.max(acc.abs()); + for &c in iter { + acc = u.mul_add(acc, c as f64 / SCALE); + worst = worst.max(acc.abs()); + worst = worst.max((u * acc).abs()); + } + } + worst + } + + #[test] + fn intermediates_stay_well_below_one() { + for table in [&SIN_P_LO[..], &SIN_P_HI[..], &COS_Q_LO[..], &COS_Q_HI[..]] { + let worst = max_intermediate(table); + assert!(worst < 0.6, "Horner intermediate reached {worst}"); + } + } + + fn saturating_horner(coeffs: &[i64; N], x: T) -> T { + let mut iter = coeffs.iter(); + let mut result = T::from_i1f63(*iter.next().unwrap()); + for &c in iter { + result = T::from_i1f63(c).saturating_add(x.saturating_mul(result)); + } + result + } + + fn check() { + for i in 0..=2000 { + let u = T::from_num(U_MAX * f64::from(i) / 2000.0); + assert_eq!(horner(&SIN_P_LO, u), saturating_horner(&SIN_P_LO, u)); + assert_eq!(horner(&SIN_P_HI, u), saturating_horner(&SIN_P_HI, u)); + assert_eq!(horner(&COS_Q_LO, u), saturating_horner(&COS_Q_LO, u)); + assert_eq!(horner(&COS_Q_HI, u), saturating_horner(&COS_Q_HI, u)); + } + } + + #[test] + fn wrapping_matches_saturating_bit_for_bit() { + check::(); + check::(); + check::(); + } +} diff --git a/tests/unit/traits.rs b/tests/unit/traits.rs index 1bb95ff..a0c556e 100644 --- a/tests/unit/traits.rs +++ b/tests/unit/traits.rs @@ -121,3 +121,106 @@ mod tests { ); } } + +/// The integer and wrapping helpers added for the fast paths. +#[cfg(test)] +#[allow(clippy::unwrap_used, reason = "test code uses unwrap for conciseness")] +mod fast_path_helpers { + use fixed::types::{I4F4, I8F8, I16F16, I32F32, I64F64}; + use fixed_analytics::CordicNumber; + + #[test] + fn div_int_matches_fixed_point_division_for_non_negative() { + for raw in [0i64, 1, 7, 12, 65_536, 100_000, i64::from(i32::MAX)] { + let x = I32F32::from_bits(raw << 20); + for k in 1..=200u32 { + assert_eq!( + x.div_int(k), + CordicNumber::div(x, I32F32::from_num(k)), + "{x} / {k}" + ); + } + } + } + + #[test] + fn div_int_truncates_toward_zero_and_saturates_on_zero() { + let x = I16F16::from_num(-1.5); + assert_eq!(x.div_int(4), I16F16::from_num(-0.375)); + assert_eq!(I16F16::from_bits(-7).div_int(2), I16F16::from_bits(-3)); + assert_eq!(I16F16::ONE.div_int(0), I16F16::MAX); + assert_eq!((-I16F16::ONE).div_int(0), I16F16::MIN); + // A divisor beyond the raw type's range exceeds |bits|. + assert_eq!(I4F4::MAX.div_int(200), I4F4::ZERO); + assert_eq!(I4F4::MIN.div_int(200), I4F4::ZERO); + } + + #[test] + fn mul_int_is_exact_and_saturates() { + assert_eq!( + I16F16::LN_2.mul_int(3), + I16F16::LN_2 + I16F16::LN_2 + I16F16::LN_2 + ); + assert_eq!(I16F16::LN_2.mul_int(-2), -(I16F16::LN_2 + I16F16::LN_2)); + assert_eq!(I16F16::LN_2.mul_int(0), I16F16::ZERO); + assert_eq!(I16F16::from_num(2).mul_int(100_000), I16F16::MAX); + assert_eq!(I16F16::from_num(-2).mul_int(100_000), I16F16::MIN); + assert_eq!(I8F8::ONE.mul_int(40_000), I8F8::MAX); + assert_eq!(I8F8::ONE.mul_int(-40_000), I8F8::MIN); + assert_eq!((-I8F8::ONE).mul_int(40_000), I8F8::MIN); + assert_eq!((-I8F8::ONE).mul_int(-40_000), I8F8::MAX); + assert_eq!(I8F8::ZERO.mul_int(40_000), I8F8::ZERO); + assert_eq!(I4F4::ONE.mul_int(200), I4F4::MAX); + } + + #[test] + fn wrapping_ops_match_saturating_ops_in_range() { + let a = I64F64::from_num(0.617); + let b = I64F64::from_num(-0.1667); + assert_eq!(a.wrapping_mul(b), a.saturating_mul(b)); + assert_eq!(a.wrapping_add(b), a.saturating_add(b)); + assert_eq!(a.wrapping_sub(b), a.saturating_sub(b)); + assert_eq!(I16F16::MAX.wrapping_add(I16F16::from_bits(1)), I16F16::MIN); + assert_eq!(I16F16::MIN.wrapping_sub(I16F16::from_bits(1)), I16F16::MAX); + } + + #[test] + fn checked_int_log2_is_floor_of_log2() { + assert_eq!(I16F16::ONE.checked_int_log2(), Some(0)); + assert_eq!(I16F16::from_num(0.5).checked_int_log2(), Some(-1)); + assert_eq!(I16F16::from_num(0.49).checked_int_log2(), Some(-2)); + assert_eq!(I16F16::from_num(6).checked_int_log2(), Some(2)); + assert_eq!(I16F16::from_num(8).checked_int_log2(), Some(3)); + assert_eq!(I16F16::from_bits(1).checked_int_log2(), Some(-16)); + assert_eq!(I16F16::ZERO.checked_int_log2(), None); + assert_eq!(I16F16::from_num(-1).checked_int_log2(), None); + assert_eq!(I64F64::MAX.checked_int_log2(), Some(62)); + } + + #[test] + fn to_i32_and_round_saturate() { + assert_eq!(I32F32::MAX.to_i32(), i32::MAX); + assert_eq!(I32F32::MIN.to_i32(), i32::MIN); + assert_eq!(I64F64::from_num(1e12).to_i32(), i32::MAX); + assert_eq!(I64F64::from_num(-1e12).to_i32(), i32::MIN); + assert_eq!(I16F16::from_num(-2.7).to_i32(), -3); + assert_eq!(I16F16::from_num(2.7).to_i32(), 2); + assert_eq!(CordicNumber::round(I16F16::MAX), I16F16::MAX); + assert_eq!(CordicNumber::round(I16F16::MIN), I16F16::MIN); + assert_eq!( + CordicNumber::round(I16F16::from_num(2.5)), + I16F16::from_num(3) + ); + assert_eq!( + CordicNumber::round(I16F16::from_num(-2.5)), + I16F16::from_num(-3) + ); + } + + #[test] + fn sqrt_round_of_negative_is_root_of_magnitude() { + assert_eq!(I16F16::from_num(-4).sqrt_round(), I16F16::from_num(2)); + assert_eq!(I64F64::from_num(-4).sqrt_round(), I64F64::from_num(2)); + assert_eq!(I16F16::ZERO.sqrt_round(), I16F16::ZERO); + } +} diff --git a/tools/accuracy-bench/baseline.json b/tools/accuracy-bench/baseline.json index 4a1181d..466ed37 100644 --- a/tools/accuracy-bench/baseline.json +++ b/tools/accuracy-bench/baseline.json @@ -1,5 +1,5 @@ { - "timestamp": 1783329055, + "timestamp": 1788385853, "results": [ { "name": "sin", @@ -95,29 +95,29 @@ "name": "asin", "i16f16": { "count": 59003, - "abs_max": 0.0001245804871132794, - "abs_mean": 0.000023135371505940558, - "abs_p50": 0.00001981861728889145, - "abs_p95": 0.000055935129377848725, - "abs_p99": 0.00007413155549662598, - "rel_max": 0.2286483921083974, - "rel_mean": 0.0001848477641879009, - "rel_p50": 0.00003599816038099953, - "rel_p95": 0.00047494272931951847, - "rel_p99": 0.0023205264290453286 + "abs_max": 0.00009864573004936261, + "abs_mean": 0.00001799490608480759, + "abs_p50": 0.000015493640990849045, + "abs_p95": 0.000043223174314654944, + "abs_p99": 0.0000570626041155875, + "rel_max": 0.19538343103455744, + "rel_mean": 0.00011305146526842695, + "rel_p50": 0.000027997800598613357, + "rel_p95": 0.0003834492838880903, + "rel_p99": 0.001275280283916785 }, "i32f32": { "count": 59003, - "abs_max": 2.64860694487723e-9, - "abs_mean": 4.734204961464557e-10, - "abs_p50": 4.015949794933249e-10, - "abs_p95": 1.1598155769121377e-9, - "abs_p99": 1.5267173636424047e-9, + "abs_max": 1.7564201204578467e-9, + "abs_mean": 3.629458328882121e-10, + "abs_p50": 3.0988378529883676e-10, + "abs_p95": 8.842215293292099e-10, + "abs_p99": 1.1424756141131809e-9, "rel_max": 6.585652811427218e-6, - "rel_mean": 4.491789537472936e-9, - "rel_p50": 7.349755952360814e-10, - "rel_p95": 8.887621226836101e-9, - "rel_p99": 4.305661252158066e-8 + "rel_mean": 3.6811745948744645e-9, + "rel_p50": 5.719933255156975e-10, + "rel_p95": 6.399402016725715e-9, + "rel_p99": 3.077119944621495e-8 }, "samples_tested": 59003 }, @@ -125,29 +125,29 @@ "name": "acos", "i16f16": { "count": 59003, - "abs_max": 0.0001345724505834589, - "abs_mean": 0.000025194457450702663, - "abs_p50": 0.000021601056192421808, - "abs_p95": 0.00006071819610875551, - "abs_p99": 0.00007890516076436427, - "rel_max": 0.0007606273037416324, - "rel_mean": 0.000023312629588548, - "rel_p50": 0.000014741930162677303, - "rel_p95": 0.00006857606436282256, - "rel_p99": 0.00018481621520396088 + "abs_max": 0.00010945006400842061, + "abs_mean": 0.000020000428168866564, + "abs_p50": 0.00001709619516243599, + "abs_p95": 0.00004841474955941116, + "abs_p99": 0.00006283908249127279, + "rel_max": 0.0005719196384412591, + "rel_mean": 0.000017890655030255327, + "rel_p50": 0.000011700957259791313, + "rel_p95": 0.000048067211811895476, + "rel_p99": 0.00015238477060454026 }, "i32f32": { "count": 59003, - "abs_max": 2.637457363618978e-9, - "abs_mean": 4.771190947426562e-10, - "abs_p50": 4.0556669134161893e-10, - "abs_p95": 1.1676429267915012e-9, - "abs_p99": 1.5423831101202268e-9, - "rel_max": 1.4374429727514298e-8, - "rel_mean": 4.513133469796149e-10, - "rel_p50": 2.711501137985831e-10, - "rel_p95": 1.3836960810872683e-9, - "rel_p99": 3.586759154363366e-9 + "abs_max": 1.8171912863351736e-9, + "abs_mean": 3.6605178957339117e-10, + "abs_p50": 3.111770841002226e-10, + "abs_p95": 8.893126235420823e-10, + "abs_p99": 1.1530281174287893e-9, + "rel_max": 1.0444955379205956e-8, + "rel_mean": 3.4956493145821594e-10, + "rel_p50": 2.1100168830853163e-10, + "rel_p95": 1.0431511349799786e-9, + "rel_p99": 2.9118122080809747e-9 }, "samples_tested": 59003 }, @@ -155,29 +155,29 @@ "name": "atan", "i16f16": { "count": 59007, - "abs_max": 0.00009635842086969104, - "abs_mean": 0.000021516147001383534, - "abs_p50": 0.000018578354218812265, - "abs_p95": 0.000051847307452890234, - "abs_p99": 0.00006198004229185372, - "rel_max": 0.00261310317711038, - "rel_mean": 0.000014991770491118831, - "rel_p50": 0.000012255592507336696, - "rel_p95": 0.00003453497462890847, - "rel_p99": 0.00004411362592735666 + "abs_max": 0.00005294208444178716, + "abs_mean": 0.000011054017025981598, + "abs_p50": 9.900200620638344e-6, + "abs_p95": 0.000026250679324713033, + "abs_p99": 0.00003439137524807734, + "rel_max": 0.0011943114886050297, + "rel_mean": 7.79437542387073e-6, + "rel_p50": 6.4947609423880385e-6, + "rel_p95": 0.000018178223294120574, + "rel_p99": 0.00002515715295884917 }, "i32f32": { "count": 59007, - "abs_max": 2.249368913354033e-9, - "abs_mean": 4.313711147886401e-10, - "abs_p50": 3.6981440132421994e-10, - "abs_p95": 1.0434733077602232e-9, - "abs_p99": 1.3546948007814308e-9, - "rel_max": 4.0305768754078694e-8, - "rel_mean": 3.00707401822319e-10, - "rel_p50": 2.4352410123890363e-10, - "rel_p95": 7.052196144414705e-10, - "rel_p99": 9.914488443352585e-10 + "abs_max": 1.4461549735500512e-9, + "abs_mean": 2.6811737611053484e-10, + "abs_p50": 2.3008150940029282e-10, + "abs_p95": 6.451517098327031e-10, + "abs_p99": 8.336746848414123e-10, + "rel_max": 4.255940723353228e-8, + "rel_mean": 1.8837865213688376e-10, + "rel_p50": 1.5136209428133091e-10, + "rel_p95": 4.382555107953913e-10, + "rel_p99": 6.327423299160185e-10 }, "samples_tested": 59007 }, @@ -305,29 +305,29 @@ "name": "asinh", "i16f16": { "count": 59007, - "abs_max": 0.009323178608780402, - "abs_mean": 0.0021284888817507635, - "abs_p50": 0.0014844659980624009, - "abs_p95": 0.006257894227960303, - "abs_p99": 0.007984492050545633, - "rel_max": 0.016260839471297594, - "rel_mean": 0.0006442542549104709, - "rel_p50": 0.0004833604162310917, - "rel_p95": 0.001753899668466157, - "rel_p99": 0.0022055481830100744 + "abs_max": 0.00018767505615269187, + "abs_mean": 0.000048199551443049315, + "abs_p50": 0.00004339764783800604, + "abs_p95": 0.00010895692598644757, + "abs_p99": 0.00013574377600322762, + "rel_max": 0.012723747609634486, + "rel_mean": 0.000024210755765449674, + "rel_p50": 0.000016135477230174978, + "rel_p95": 0.000050803824758986083, + "rel_p99": 0.00010673290949153741 }, "i32f32": { "count": 59007, - "abs_max": 1.765475614590173e-7, - "abs_mean": 3.378897481685445e-8, - "abs_p50": 2.2932562071531493e-8, - "abs_p95": 1.0207803402551008e-7, - "abs_p99": 1.309457595688457e-7, - "rel_max": 4.6586130728043937e-7, - "rel_mean": 1.0284023966828364e-8, - "rel_p50": 7.592837972471504e-9, - "rel_p95": 2.849768247469795e-8, - "rel_p99": 3.622677254418828e-8 + "abs_max": 4.528247998791812e-9, + "abs_mean": 1.412659747773413e-9, + "abs_p50": 1.3856520375554737e-9, + "abs_p95": 2.770816465158532e-9, + "abs_p99": 3.329308828625699e-9, + "rel_max": 3.68994826076661e-7, + "rel_mean": 6.273296847090133e-10, + "rel_p50": 5.007130347601334e-10, + "rel_p95": 1.1890082070195224e-9, + "rel_p99": 2.391133491764771e-9 }, "samples_tested": 59007 }, @@ -335,29 +335,29 @@ "name": "acosh", "i16f16": { "count": 59001, - "abs_max": 0.009688040314605129, - "abs_mean": 0.0022273528018676597, - "abs_p50": 0.0015914165377841627, - "abs_p95": 0.006447086883274444, - "abs_p99": 0.008295719356268272, - "rel_max": 0.002626733490764648, - "rel_mean": 0.0006736969672385671, - "rel_p50": 0.0005214367925851396, - "rel_p95": 0.0018006676087912823, - "rel_p99": 0.0022695168421162757 + "abs_max": 0.00018582392440702478, + "abs_mean": 0.000045107839795604246, + "abs_p50": 0.00004015930818712654, + "abs_p95": 0.00010417120364358823, + "abs_p99": 0.00013155016010113485, + "rel_max": 0.0006137833808951901, + "rel_mean": 0.000018746843532914047, + "rel_p50": 0.000014809569969610293, + "rel_p95": 0.00004654927607609366, + "rel_p99": 0.00008639361951765569 }, "i32f32": { "count": 59001, - "abs_max": 1.7263420248880834e-7, - "abs_mean": 3.481698610113588e-8, - "abs_p50": 2.4109737317701274e-8, - "abs_p95": 1.0343075285135228e-7, - "abs_p99": 1.3352229855101427e-7, - "rel_max": 4.68161453820617e-8, - "rel_mean": 1.0532948305768858e-8, - "rel_p50": 7.964701008948365e-9, - "rel_p95": 2.8769664915903364e-8, - "rel_p99": 3.6546740078326744e-8 + "abs_max": 4.775437378867764e-9, + "abs_mean": 1.3471755089957144e-9, + "abs_p50": 1.30687816124464e-9, + "abs_p95": 2.6884223736090007e-9, + "abs_p99": 3.2257703175275765e-9, + "rel_max": 1.5236174923360043e-8, + "rel_mean": 5.290218757158357e-10, + "rel_p50": 4.749859389075313e-10, + "rel_p95": 1.1120435645793785e-9, + "rel_p99": 1.8034071434472156e-9 }, "samples_tested": 59001 }, @@ -365,29 +365,29 @@ "name": "atanh", "i16f16": { "count": 59001, - "abs_max": 0.000996936539403137, - "abs_mean": 0.000048225767174518445, - "abs_p50": 0.00003234791334638665, - "abs_p95": 0.00013384945319838693, - "abs_p99": 0.00036007364172752077, - "rel_max": 0.3650807884541502, - "rel_mean": 0.0003009640696123398, - "rel_p50": 0.0000590390888308655, - "rel_p95": 0.0006253162143217303, - "rel_p99": 0.003043533591188758 + "abs_max": 0.001042712906590637, + "abs_mean": 0.000036939965338882025, + "abs_p50": 0.000020152732652878314, + "abs_p95": 0.00012821044868305265, + "abs_p99": 0.0003632297675979501, + "rel_max": 0.39531509631821815, + "rel_mean": 0.00022127427735330204, + "rel_p50": 0.00003985094006796104, + "rel_p95": 0.0003334768129343867, + "rel_p99": 0.0016280093761588796 }, "i32f32": { "count": 59001, - "abs_max": 1.7351338144067086e-8, - "abs_mean": 1.0091975489373658e-9, - "abs_p50": 7.656760780960781e-10, - "abs_p95": 2.5722681762374577e-9, - "abs_p99": 5.487656551395048e-9, - "rel_max": 8.614530652505506e-6, - "rel_mean": 6.677347018196235e-9, - "rel_p50": 1.322534893292658e-9, - "rel_p95": 1.443044339663817e-8, - "rel_p99": 7.927532805879432e-8 + "abs_max": 1.6187184925797737e-8, + "abs_mean": 6.261754052554837e-10, + "abs_p50": 3.836353457131736e-10, + "abs_p95": 1.9506158910331806e-9, + "abs_p99": 5.410538683747745e-9, + "rel_max": 6.587852476634234e-6, + "rel_mean": 3.633106045916069e-9, + "rel_p50": 7.203825977353904e-10, + "rel_p95": 6.242463909235762e-9, + "rel_p99": 3.2775443321514364e-8 }, "samples_tested": 59001 }, @@ -395,29 +395,29 @@ "name": "acoth", "i16f16": { "count": 59001, - "abs_max": 0.00063469181628939, - "abs_mean": 0.00003892003359818683, - "abs_p50": 0.00003349354136025079, - "abs_p95": 0.00009199005491542477, - "abs_p99": 0.00011657816953467709, - "rel_max": 0.01284496346263612, - "rel_mean": 0.0020995023707175223, - "rel_p50": 0.0013292638519312421, - "rel_p95": 0.006665326707366214, - "rel_p99": 0.009055939297598636 + "abs_max": 0.000644926102925325, + "abs_mean": 0.000022406057733804776, + "abs_p50": 0.000019721135377757937, + "abs_p95": 0.00005041140596961924, + "abs_p99": 0.00006179515435000837, + "rel_max": 0.006806962761098952, + "rel_mean": 0.0012058340747461056, + "rel_p50": 0.0006845795226896409, + "rel_p95": 0.004087593711076304, + "rel_p99": 0.005224262624724212 }, "i32f32": { "count": 59001, - "abs_max": 1.008287364712146e-8, - "abs_mean": 8.210808801982123e-10, - "abs_p50": 6.99425746486515e-10, - "abs_p95": 2.000697003901042e-9, - "abs_p99": 2.6134055180343507e-9, - "rel_max": 4.2790418218212674e-7, - "rel_mean": 4.255886758338163e-8, - "rel_p50": 2.617976365074166e-8, - "rel_p95": 1.3923560561915121e-7, - "rel_p99": 2.0740532514525208e-7 + "abs_max": 1.0548534934429199e-8, + "abs_mean": 3.868739513967378e-10, + "abs_p50": 3.349050689549493e-10, + "abs_p95": 9.021436098155533e-10, + "abs_p99": 1.1805792975161378e-9, + "rel_max": 1.501734919624401e-7, + "rel_mean": 1.9445563374156688e-8, + "rel_p50": 1.2266884421399114e-8, + "rel_p95": 6.27902243930849e-8, + "rel_p99": 8.86609767234068e-8 }, "samples_tested": 59001 }, @@ -455,29 +455,29 @@ "name": "ln", "i16f16": { "count": 59003, - "abs_max": 0.007120513357136815, - "abs_mean": 0.00006108372519995839, - "abs_p50": 0.0000507042797313062, - "abs_p95": 0.0001528340420877683, - "abs_p99": 0.00020458018328906036, - "rel_max": 0.01244179019784243, - "rel_mean": 0.000013511077843052894, - "rel_p50": 8.759334514528923e-6, - "rel_p95": 0.00002967197208423552, - "rel_p99": 0.00005295682216701669 + "abs_max": 0.007059478200886815, + "abs_mean": 0.000048678384609382866, + "abs_p50": 0.00004412568314293708, + "abs_p95": 0.00010925010635975951, + "abs_p99": 0.00013469276996680435, + "rel_max": 0.08319029452329336, + "rel_mean": 0.000011832575993244523, + "rel_p50": 7.614814836016055e-6, + "rel_p95": 0.000020924425418536253, + "rel_p99": 0.0000369448686848002 }, "i32f32": { "count": 59003, - "abs_max": 6.868512514301983e-8, - "abs_mean": 2.2163491276603115e-9, - "abs_p50": 1.9963826147773034e-9, - "abs_p95": 5.00238517275875e-9, - "abs_p99": 6.2628524588603796e-9, - "rel_max": 7.379229166444631e-7, - "rel_mean": 4.50470369587187e-10, - "rel_p50": 3.4815093893310755e-10, - "rel_p95": 9.17042603709098e-10, - "rel_p99": 1.4470851177100795e-9 + "abs_max": 6.821946385571209e-8, + "abs_mean": 2.205632261181594e-9, + "abs_p50": 2.2086270590193635e-9, + "abs_p95": 3.6060017194472493e-9, + "abs_p99": 4.177266532678914e-9, + "rel_max": 1.5151121336623418e-6, + "rel_mean": 4.3492454435577516e-10, + "rel_p50": 3.7791581174422563e-10, + "rel_p95": 6.473037250784497e-10, + "rel_p99": 9.450187564328677e-10 }, "samples_tested": 59003 }, @@ -485,29 +485,29 @@ "name": "log2", "i16f16": { "count": 59003, - "abs_max": 0.0006598014362131366, - "abs_mean": 0.00008558332538114927, - "abs_p50": 0.00007083049510292483, - "abs_p95": 0.00021631260077370484, - "abs_p99": 0.00028929750223483097, + "abs_max": 0.0007055778034006366, + "abs_mean": 0.00006387910309026122, + "abs_p50": 0.00005647667593500216, + "abs_p95": 0.00014796261375948916, + "abs_p99": 0.0001858945163277781, "rel_max": 0.014636494762760797, - "rel_mean": 0.000013259536025697493, - "rel_p50": 8.478787777090254e-6, - "rel_p95": 0.000029151110902145505, - "rel_p99": 0.00005287799348650687 + "rel_mean": 9.980888854233613e-6, + "rel_p50": 6.720821630152601e-6, + "rel_p95": 0.000020138496604533648, + "rel_p99": 0.00003782534602321454 }, "i32f32": { "count": 59003, - "abs_max": 1.21446355194621e-8, - "abs_mean": 2.203859919803808e-9, - "abs_p50": 1.8789663158713665e-9, - "abs_p95": 5.352216447818137e-9, - "abs_p99": 6.9875891739457074e-9, - "rel_max": 5.933103149548927e-7, - "rel_mean": 3.4640207161669755e-10, - "rel_p50": 2.241559775657736e-10, - "rel_p95": 7.212274150853486e-10, - "rel_p99": 1.319858054780029e-9 + "abs_max": 6.775742633635673e-9, + "abs_mean": 1.2482876363506192e-9, + "abs_p50": 1.1024754442701123e-9, + "abs_p95": 2.9128237599707063e-9, + "abs_p99": 3.64468100144677e-9, + "rel_max": 3.207643232394731e-7, + "rel_mean": 1.9224994832014074e-10, + "rel_p50": 1.3147901067393455e-10, + "rel_p95": 3.9040981853032305e-10, + "rel_p99": 7.121447298153111e-10 }, "samples_tested": 59003 }, @@ -516,28 +516,28 @@ "i16f16": { "count": 59003, "abs_max": 0.0001983642578125, - "abs_mean": 0.00002802985157528529, - "abs_p50": 0.000023431607028001622, - "abs_p95": 0.00007010842908883319, - "abs_p99": 0.000092947326506998, + "abs_mean": 0.000024233678656285957, + "abs_p50": 0.00002262463298530193, + "abs_p95": 0.00005244047921992845, + "abs_p99": 0.0000638226606906045, "rel_max": 0.020974729725972013, - "rel_mean": 0.000014423993800471001, - "rel_p50": 9.28207260892599e-6, - "rel_p95": 0.000031374333977432835, - "rel_p99": 0.00005676337528031852 + "rel_mean": 0.000012491983123135089, + "rel_p50": 8.958415045300318e-6, + "rel_p95": 0.000023629069645499613, + "rel_p99": 0.00004474934477675163 }, "i32f32": { "count": 59003, - "abs_max": 4.247256324418913e-9, - "abs_mean": 9.19977012467165e-10, - "abs_p50": 8.189431355276611e-10, - "abs_p95": 2.117753084007745e-9, - "abs_p99": 2.6576736367189824e-9, - "rel_max": 5.877489078123925e-7, - "rel_mean": 4.485086856900398e-10, - "rel_p50": 3.2691886929886964e-10, - "rel_p95": 9.069053951197586e-10, - "rel_p99": 1.5008492405980964e-9 + "abs_max": 2.4195725423226122e-9, + "abs_mean": 9.018912722648107e-10, + "abs_p50": 9.034386572182029e-10, + "abs_p95": 1.5063625902200783e-9, + "abs_p99": 1.7491261772306643e-9, + "rel_max": 4.368525433076608e-7, + "rel_mean": 4.0554272131845393e-10, + "rel_p50": 3.5657916451437263e-10, + "rel_p95": 6.390202251478334e-10, + "rel_p99": 1.006565285033717e-9 }, "samples_tested": 59003 }, @@ -575,29 +575,29 @@ "name": "sqrt", "i16f16": { "count": 59003, - "abs_max": 0.000018870804664450347, - "abs_mean": 7.633058381899692e-6, - "abs_p50": 7.615373824876315e-6, - "abs_p95": 0.000014475517588152798, - "abs_p99": 0.000015094551585548288, - "rel_max": 0.00004366209963999045, - "rel_mean": 1.768162523538784e-7, - "rel_p50": 1.1613413047229915e-7, - "rel_p95": 4.73581011416962e-7, - "rel_p99": 1.3197749510934042e-6 + "abs_max": 8.484410951192789e-6, + "abs_mean": 3.82567097197237e-6, + "abs_p50": 3.832319862340228e-6, + "abs_p95": 7.246182150311142e-6, + "abs_p99": 7.550143081402894e-6, + "rel_max": 0.000024877819481652173, + "rel_mean": 8.883583091079036e-8, + "rel_p50": 5.802056761914492e-8, + "rel_p95": 2.4228606832629987e-7, + "rel_p99": 6.764313336895449e-7 }, "i32f32": { "count": 59003, - "abs_max": 3.707159579313668e-10, - "abs_mean": 1.165586754832912e-10, - "abs_p50": 1.1648637610051082e-10, - "abs_p95": 2.211919536421192e-10, - "abs_p99": 2.304574309164309e-10, - "rel_max": 1.3218584815523091e-9, - "rel_mean": 2.6958073035473984e-12, - "rel_p50": 1.776147695366651e-12, - "rel_p95": 7.161895644954702e-12, - "rel_p99": 2.0013121754077194e-11 + "abs_max": 1.7630552573422165e-10, + "abs_mean": 5.8168412384479286e-11, + "abs_p50": 5.800870894745458e-11, + "abs_p95": 1.1058887139370199e-10, + "abs_p99": 1.1536727129168867e-10, + "rel_max": 6.254206431271312e-10, + "rel_mean": 1.370321034827564e-12, + "rel_p50": 8.848455916612933e-13, + "rel_p95": 3.620677038772172e-12, + "rel_p99": 1.0162809373904082e-11 }, "samples_tested": 59003 }