From a551395f0887f4a77e33038527ebf25b33b92648 Mon Sep 17 00:00:00 2001 From: GeEom Date: Thu, 3 Sep 2026 09:33:30 +0100 Subject: [PATCH 1/4] Replace fixed-point divisions where results are unchanged MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A 128-bit fixed-point division costs about twelve multiplications, and forming a divisor like 182 in T overflows types with few integer bits. CordicNumber gains div_int and mul_int, which work on the raw representation, plus wrapping_mul/add/sub and checked_int_log2; round and to_i32 now saturate instead of wrapping at MAX. exp and sinh_cosh divide by their Taylor divisors with div_int. sin_cos and exp estimate their reductions with stored reciprocals when the integer bits do not outnumber the fractional bits, and correct the estimate once in either direction. ln normalises from the leading bit instead of a shift loop. Horner evaluation wraps, since every intermediate stays below 1 for the tables in use. Results are bit-identical on every layout with at least 8 integer bits except where the old code was wrong: exp above 2^31·ln 2 on 64-bit and wider types wrapped in to_i32, and angles within π of MAX saturated the reduction to sin_cos(0). I4F12, I4F60 and I3F125 gain working exp and sinh_cosh. I64F64: sin_cos 52 -> 26 ns, exp 159 -> 41 ns, sinh_cosh 175 -> 54 ns, and ln normalisation is O(1). Adds I64F64 benchmarks, instantiates the no-panic binary on three layouts, and raises the no-panic floor to 0.1.37. --- Cargo.toml | 2 +- benches/benchmarks.rs | 115 ++++++++++++++----------- src/bin/verify_no_panic.rs | 49 ++++++++--- src/ops/circular.rs | 24 +++++- src/ops/exponential.rs | 126 +++++++++++++++------------ src/ops/hyperbolic.rs | 58 +++++-------- src/tables/chebyshev.rs | 7 +- src/traits.rs | 95 +++++++++++++++++++-- tests/unit/mod.rs | 1 + tests/unit/ops/circular.rs | 73 ++++++++++++++++ tests/unit/ops/exponential.rs | 151 +++++++++++++++++++++++++++++++++ tests/unit/ops/hyperbolic.rs | 32 +++++++ tests/unit/support.rs | 29 +++++++ tests/unit/tables/chebyshev.rs | 67 +++++++++++++++ tests/unit/traits.rs | 96 +++++++++++++++++++++ 15 files changed, 755 insertions(+), 170 deletions(-) create mode 100644 tests/unit/support.rs diff --git a/Cargo.toml b/Cargo.toml index a1044c0..a46f6bd 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -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/benches/benchmarks.rs b/benches/benchmarks.rs index a8786f5..9bf3dec 100644 --- a/benches/benchmarks.rs +++ b/benches/benchmarks.rs @@ -1,68 +1,81 @@ -//! 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, 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.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/src/bin/verify_no_panic.rs b/src/bin/verify_no_panic.rs index 6d1c735..8e97e4e 100644 --- a/src/bin/verify_no_panic.rs +++ b/src/bin/verify_no_panic.rs @@ -1,13 +1,16 @@ -//! 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; @@ -16,11 +19,7 @@ use fixed_analytics::{ 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,35 @@ 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)); // 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/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..426d175 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,69 @@ 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(); - - // 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 { + // 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(); + + // 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 and 1.4e-12 + // at degree 12: below half an ulp up to 18 and 38 fractional bits. let mut p = one; 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)); + 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)); } - // 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(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" @@ -165,27 +180,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)) diff --git a/src/ops/hyperbolic.rs b/src/ops/hyperbolic.rs index fe175e1..19ebc27 100644 --- a/src/ops/hyperbolic.rs +++ b/src/ops/hyperbolic.rs @@ -64,44 +64,28 @@ 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 + // and 3.7e-12 / 2.9e-13 at 13 / 14: below half an ulp up to 22 and 37 + // 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() >= 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) 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..fff2ad0 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -90,19 +90,44 @@ 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; /// 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; } @@ -240,6 +265,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 +333,11 @@ macro_rules! impl_cordic_generic { } } + #[inline] + fn checked_int_log2(self) -> Option { + Fixed::checked_int_log2(self) + } + #[inline] fn from_num(n: N) -> Self { Self::from_num(n) @@ -272,16 +355,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::() } } }; 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/circular.rs b/tests/unit/ops/circular.rs index 9de8625..6d8c1eb 100644 --- a/tests/unit/ops/circular.rs +++ b/tests/unit/ops/circular.rs @@ -516,3 +516,76 @@ 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, 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() + ); + } + } +} diff --git a/tests/unit/ops/exponential.rs b/tests/unit/ops/exponential.rs index b62a74a..b9e013c 100644 --- a/tests/unit/ops/exponential.rs +++ b/tests/unit/ops/exponential.rs @@ -424,3 +424,154 @@ 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, pow2}; + + #[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(); + // Degree-12 Taylor: truncation ≈ r¹³/13! ≈ 1.4e-12 of the result. + let tol = 8.0f64.mul_add(2f64.powi(-64), 3e-12 * 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() < 3e-12, + "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 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..561717c 100644 --- a/tests/unit/ops/hyperbolic.rs +++ b/tests/unit/ops/hyperbolic.rs @@ -487,3 +487,35 @@ 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 fixed::types::{I4F12, I4F60}; + use fixed_analytics::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() < 1e-11, + "I4F60 sinh({v60}) = {s60}" + ); + assert!( + (c60 - v60.cosh()).abs() < 1e-11, + "I4F60 cosh({v60}) = {c60}" + ); + } + } +} 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..5948165 100644 --- a/tests/unit/traits.rs +++ b/tests/unit/traits.rs @@ -121,3 +121,99 @@ 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) + ); + } +} From 917ec322fc707cdbdc4f65809b05ffb6a0365cb5 Mon Sep 17 00:00:00 2001 From: GeEom Date: Thu, 3 Sep 2026 09:33:34 +0100 Subject: [PATCH 2/4] Add pow MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit pow(base, exponent) = exp(exponent · ln(base)), fallible on a negative base or a zero base with a negative exponent; 0^0 is 1 and exp's saturation applies. ln is split into its domain check and ln_positive so pow can skip the check. The logarithm's rounding error is multiplied by the exponent, so the relative error grows with |exponent|: about 1e-4 at I16F16 for exponent 1.5. 219 ns on I64F64, 29 ns on I16F16. --- README.md | 2 +- benches/benchmarks.rs | 5 +++- src/bin/verify_no_panic.rs | 3 ++- src/lib.rs | 4 +-- src/ops/exponential.rs | 43 +++++++++++++++++++++++++++----- src/ops/mod.rs | 4 +-- tests/unit/ops/exponential.rs | 47 ++++++++++++++++++++++++++++++++++- 7 files changed, 94 insertions(+), 14 deletions(-) diff --git a/README.md b/README.md index 62f1357..a2002ed 100644 --- a/README.md +++ b/README.md @@ -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. diff --git a/benches/benchmarks.rs b/benches/benchmarks.rs index 9bf3dec..399fa3a 100644 --- a/benches/benchmarks.rs +++ b/benches/benchmarks.rs @@ -8,7 +8,7 @@ use criterion::{Criterion, criterion_group, criterion_main}; use fixed::types::{I16F16, I64F64}; use fixed_analytics::{ CordicNumber, 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, + log2, log10, pow, pow2, sin, sin_cos, sinh, sinh_cosh, sqrt, tan, tanh, }; fn bench_type(c: &mut Criterion, name: &str) { @@ -58,6 +58,9 @@ fn bench_type(c: &mut Criterion, name: &str) { 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(); } { diff --git a/src/bin/verify_no_panic.rs b/src/bin/verify_no_panic.rs index 8e97e4e..47d38c9 100644 --- a/src/bin/verify_no_panic.rs +++ b/src/bin/verify_no_panic.rs @@ -16,7 +16,7 @@ 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 exercise(x: T, y: T, two: T) { @@ -47,6 +47,7 @@ fn exercise(x: T, y: T, two: T) { let _ = std::hint::black_box(atanh(x)); let _ = std::hint::black_box(coth(x)); let _ = std::hint::black_box(acoth(two)); + let _ = std::hint::black_box(pow(two, x)); // Type-safe wrapper functions if let Some(nn) = NonNegative::new(x) { 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/exponential.rs b/src/ops/exponential.rs index 426d175..750765e 100644 --- a/src/ops/exponential.rs +++ b/src/ops/exponential.rs @@ -156,6 +156,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 @@ -163,16 +190,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: @@ -215,7 +246,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/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/tests/unit/ops/exponential.rs b/tests/unit/ops/exponential.rs index b9e013c..f24ada6 100644 --- a/tests/unit/ops/exponential.rs +++ b/tests/unit/ops/exponential.rs @@ -431,7 +431,7 @@ mod tests { mod reduction { use crate::unit::support::Lcg; use fixed::types::{I4F60, I16F16, I24F8, I32F32, I64F64}; - use fixed_analytics::{CordicNumber, exp, ln, log2, pow2}; + use fixed_analytics::{CordicNumber, exp, ln, log2, pow, pow2, sqrt}; #[test] fn exp_saturates_at_the_extremes_of_wide_types() { @@ -561,6 +561,51 @@ mod reduction { 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() < 2.0, "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), 3e-12 * 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-13, "2.5^0.5 = {root}, want {want}"); + let inv: f64 = pow(x, -I64F64::ONE).unwrap().to_num(); + assert!((inv - 0.4).abs() < 1e-13, "2.5^-1 = {inv}"); + } + #[test] fn ln_tracks_f64_at_i64f64() { let mut rng = Lcg(0x1A); From 75b64af1cb4c4f10edf9799ecb0dfd2e3d7290fb Mon Sep 17 00:00:00 2001 From: GeEom Date: Thu, 3 Sep 2026 09:33:38 +0100 Subject: [PATCH 3/4] Make sqrt exact, raise wide-type precision, bump to 3.0.0 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit sqrt refined a seed with an absolute 2^-32 convergence epsilon, so on I64F64 it was exact above 100 and returned relative errors above 1 below 10^-13. It now returns the representable value nearest the true root at every magnitude: the integer root of the raw bits, plus one Newton step and a 256-bit remainder test on 128-bit types. Verified against the fixed crate's rounded-down root on 14 layouts. I64F64 sqrt 114 -> 42 ns above 1 and 180 -> 24 ns below. CORDIC vectoring stops after ⌈frac/3⌉ shift stages and closes the residual with one division; the skipped stages only added table angles of exactly 2^-i. atan2 212 -> 66 ns and ln 240 -> 127 ns on I64F64, and the hyperbolic kernel delivers 64 bits instead of 54. asinh and acosh use the logarithmic form above 1 and 1.5, where the atanh argument reduction tripled the error at every step. exp and sinh_cosh gain a third series tier for 40 or more fractional bits (degree 19, and 21/22), taking their I64F64 relative error from 1e-13 to 2e-18; the previous comment understated the degree-12 truncation by 400x. Accuracy gate against the 2.1.0 baseline: 30 improve, 10 unchanged, none regress. asinh I16F16 mean 6.44e-4 -> 2.42e-5, sqrt I32F32 2.70e-12 -> 1.37e-12, atan I16F16 1.50e-5 -> 7.79e-6. CordicNumber gained required methods and the kernels now return the residual vector, so this is a major version. --- Cargo.toml | 2 +- README.md | 26 +- src/bounded.rs | 16 ++ src/kernel/cordic.rs | 128 ++++----- src/ops/algebraic.rs | 87 +----- src/ops/exponential.rs | 14 +- src/ops/hyperbolic.rs | 91 ++++-- src/traits.rs | 182 +++++++++++- tests/unit/kernel/cordic.rs | 51 ++++ tests/unit/ops/algebraic.rs | 214 ++++++++++++++ tests/unit/ops/circular.rs | 28 +- tests/unit/ops/exponential.rs | 14 +- tests/unit/ops/hyperbolic.rs | 86 +++++- tests/unit/traits.rs | 7 + tools/accuracy-bench/baseline.json | 434 ++++++++++++++--------------- 15 files changed, 959 insertions(+), 421 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index a46f6bd..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"] diff --git a/README.md b/README.md index a2002ed..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 @@ -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/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/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/exponential.rs b/src/ops/exponential.rs index 750765e..f90ff59 100644 --- a/src/ops/exponential.rs +++ b/src/ops/exponential.rs @@ -116,9 +116,19 @@ pub fn exp(x: T) -> T { } // 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 and 1.4e-12 - // at degree 12: below half an ulp up to 18 and 38 fractional bits. + // 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 { p = one.saturating_add(r.div_int(12).saturating_mul(p)); p = one.saturating_add(r.div_int(11).saturating_mul(p)); diff --git a/src/ops/hyperbolic.rs b/src/ops/hyperbolic.rs index 19ebc27..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. @@ -65,11 +66,21 @@ pub fn sinh_cosh(x: T) -> (T, T) { let u = reduced.saturating_mul(reduced); // 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 - // and 3.7e-12 / 2.9e-13 at 13 / 14: below half an ulp up to 22 and 37 - // fractional bits respectively. + // 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 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)); @@ -173,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`. @@ -195,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/traits.rs b/src/traits.rs index fff2ad0..8a465c3 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -116,6 +116,10 @@ pub trait CordicNumber: /// 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. @@ -132,6 +136,126 @@ pub trait CordicNumber: 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 // ============================================================================= @@ -147,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) @@ -338,6 +464,11 @@ macro_rules! impl_cordic_generic { 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) @@ -377,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/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 6d8c1eb..4338348 100644 --- a/tests/unit/ops/circular.rs +++ b/tests/unit/ops/circular.rs @@ -523,7 +523,7 @@ mod tests { mod reduction_and_accuracy { use crate::unit::support::Lcg; use fixed::types::{I24F8, I32F32, I64F64}; - use fixed_analytics::{CordicNumber, sin_cos}; + 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 { @@ -588,4 +588,30 @@ mod reduction_and_accuracy { ); } } + + #[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 f24ada6..d220795 100644 --- a/tests/unit/ops/exponential.rs +++ b/tests/unit/ops/exponential.rs @@ -482,8 +482,8 @@ mod reduction { let x = I64F64::from_num(v); let got: f64 = exp(x).to_num(); let want = x.to_num::().exp(); - // Degree-12 Taylor: truncation ≈ r¹³/13! ≈ 1.4e-12 of the result. - let tol = 8.0f64.mul_add(2f64.powi(-64), 3e-12 * want); + // 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}"); } } @@ -497,7 +497,7 @@ mod reduction { let got: f64 = exp(x).to_num(); let want = x.to_num::().exp(); assert!( - ((got - want) / want).abs() < 3e-12, + ((got - want) / want).abs() < 5e-15, "I4F60 exp({v}) = {got}, want {want}" ); } @@ -579,7 +579,7 @@ mod reduction { let got: f64 = pow(I16F16::from_num(2), I16F16::from_num(10)) .unwrap() .to_num(); - assert!((got - 1024.0).abs() < 2.0, "I16F16 2^10 = {got}"); + assert!((got - 1024.0).abs() < 0.5, "I16F16 2^10 = {got}"); } #[test] @@ -592,7 +592,7 @@ mod reduction { 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), 3e-12 * want); + let tol = 8.0f64.mul_add(2f64.powi(-64), 1e-14 * want); assert!( (got - want).abs() < tol, "pow({b}, {e}) = {got}, want {want}" @@ -601,9 +601,9 @@ mod reduction { 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-13, "2.5^0.5 = {root}, want {want}"); + 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-13, "2.5^-1 = {inv}"); + assert!((inv - 0.4).abs() < 1e-15, "2.5^-1 = {inv}"); } #[test] diff --git a/tests/unit/ops/hyperbolic.rs b/tests/unit/ops/hyperbolic.rs index 561717c..0dc0f65 100644 --- a/tests/unit/ops/hyperbolic.rs +++ b/tests/unit/ops/hyperbolic.rs @@ -492,8 +492,9 @@ mod tests { #[cfg(test)] #[allow(clippy::unwrap_used, reason = "test code uses unwrap for conciseness")] mod wide_and_narrow { - use fixed::types::{I4F12, I4F60}; - use fixed_analytics::sinh_cosh; + 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() { @@ -509,13 +510,90 @@ mod wide_and_narrow { 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() < 1e-11, + (s60 - v60.sinh()).abs() < 4e-15, "I4F60 sinh({v60}) = {s60}" ); assert!( - (c60 - v60.cosh()).abs() < 1e-11, + (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/traits.rs b/tests/unit/traits.rs index 5948165..a0c556e 100644 --- a/tests/unit/traits.rs +++ b/tests/unit/traits.rs @@ -216,4 +216,11 @@ mod fast_path_helpers { 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 } From c4ddb488ed370b15aa9cb1d9a8699fcb02848641 Mon Sep 17 00:00:00 2001 From: GeEom Date: Thu, 3 Sep 2026 09:48:57 +0100 Subject: [PATCH 4/4] Exclude untestable targets from coverage Codecov's ignore pattern did not match the bench file, which sits directly under benches/, and the no-panic binary is compiled only under its feature, so neither can be covered by tests; both are excluded by folder. --- codecov.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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"