v0.5.0: Tune estimated thresholds via accuracy search

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "fixed_analytics"
-version = "0.4.0"
+version = "0.5.0"
 edition = "2024"
 rust-version = "1.88"
 authors = ["David Gathercole"]

README.md

@@ -30,33 +30,53 @@ Requires Rust 1.88 or later.
 
 ```toml
 [dependencies]
-fixed_analytics = "0.4"
+fixed_analytics = "0.5"
 ```
 
 For `no_std` environments:
 
 ```toml
 [dependencies]
-fixed_analytics = { version = "0.4", default-features = false }
+fixed_analytics = { version = "0.5", default-features = false }
 ```
 
 ## Available Functions
 
+### Function Categories
+
+**Total functions** return `T` directly and handle all inputs, possibly with saturation.
+**Fallible functions** return `Result<T, Error>` and fail on domain violations.
+
 | Category | Total Functions | Fallible Functions |
 |----------|-----------------|-------------------|
 | Trigonometric | `sin`, `cos`, `tan`, `sin_cos`, `atan`, `atan2` | `asin`, `acos` |
-| Hyperbolic | `sinh`, `cosh`, `tanh`, `sinh_cosh`, `asinh` | `acosh`, `atanh`, `acoth`, `coth` |
+| Hyperbolic | `sinh`, `cosh`, `tanh`, `sinh_cosh`\*, `asinh` | `acosh`, `atanh`, `acoth`, `coth` |
 | Exponential | `exp`, `pow2` | `ln`, `log2`, `log10` |
 | Algebraic | — | `sqrt` |
 
-Fallible functions return `Result<T, Error>` and fail on domain violations.
-Total functions return `T` directly and handle all inputs.
 Functions are calculated via CORDIC, Newton-Raphson, and Taylor series techniques.
 
+### Saturation Behavior
+
+The following total functions use **saturation semantics**: rather than returning an error or
+panicking on overflow, they clamp to the representable range.
+
+| Function | I16F16 Threshold | I32F32 Threshold | Result |
+|----------|------------------|------------------|--------|
+| `exp` | x ≥ 22.2 | x ≥ 44.4 | `T::MAX` |
+| `exp` | x ≤ -9.2 | x ≤ -16.2 | Zero |
+| `pow2` | x ≥ 15.0 | x ≥ 31.0 | `T::MAX` |
+| `pow2` | x ≤ -13.2 | x ≤ -23.3 | Zero |
+| `sinh` | \|x\| ≥ 11.1 | \|x\| ≥ 22.2 | `T::MAX` or `T::MIN` |
+| `cosh` | \|x\| ≥ 11.1 | \|x\| ≥ 22.2 | `T::MAX` |
+| `tan` | \|x - pole\| < 8e-5 | \|x - pole\| < 5e-10 | `T::MAX` or `T::MIN` |
+
+Where for `tan`, "pole" refers to ±π/2, ±3π/2, ±5π/2, ...
+
 <!-- ACCURACY_START -->
 ### Accuracy
 
-Relative error statistics measured against MPFR reference implementations.
+Relative error statistics measured against MPFR reference implementations. The file tools/accuracy-bench/baseline.json contains further measurements.
 
 | Function | I16F16 Mean | I16F16 Median | I16F16 P95 | I32F32 Mean | I32F32 Median | I32F32 P95 |
 |----------|-------------|---------------|------------|-------------|---------------|------------|
@@ -66,18 +86,18 @@ Relative error statistics measured against MPFR reference implementations.
 | asin | 2.87e-4 | 5.93e-5 | 6.46e-4 | 5.34e-9 | 8.82e-10 | 1.03e-8 |
 | acos | 3.61e-5 | 2.18e-5 | 1.14e-4 | 5.37e-10 | 3.19e-10 | 1.71e-9 |
 | atan | 2.71e-5 | 2.21e-5 | 6.29e-5 | 3.69e-10 | 2.92e-10 | 8.74e-10 |
-| sinh | 1.85e-4 | 1.39e-4 | 5.05e-4 | 1.05e-8 | 2.37e-9 | 9.47e-9 |
-| cosh | 1.73e-4 | 1.23e-4 | 5.00e-4 | 1.03e-8 | 2.07e-9 | 9.16e-9 |
-| tanh | 2.26e-5 | 1.38e-5 | 6.13e-5 | 1.67e-9 | 1.31e-10 | 1.22e-9 |
-| coth | 1.74e-5 | 4.86e-6 | 7.23e-5 | 4.34e-10 | 1.39e-10 | 1.31e-9 |
+| sinh | 1.82e-4 | 1.34e-4 | 5.03e-4 | 1.05e-8 | 2.35e-9 | 9.30e-9 |
+| cosh | 1.73e-4 | 1.23e-4 | 5.00e-4 | 1.02e-8 | 2.07e-9 | 9.16e-9 |
+| tanh | 2.08e-5 | 1.38e-5 | 5.89e-5 | 1.64e-9 | 1.31e-10 | 1.23e-9 |
+| coth | 1.20e-5 | 4.83e-6 | 5.57e-5 | 4.00e-10 | 1.39e-10 | 1.30e-9 |
 | 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 |
-| exp | 1.14e-2 | 7.72e-5 | 7.87e-2 | 1.91e-7 | 2.35e-9 | 1.30e-6 |
+| exp | 1.14e-2 | 6.67e-5 | 7.87e-2 | 1.91e-7 | 2.24e-9 | 1.30e-6 |
 | 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 |
-| pow2 | 7.30e-4 | 5.68e-5 | 4.70e-3 | 1.15e-8 | 1.16e-9 | 7.38e-8 |
+| pow2 | 7.28e-4 | 4.74e-5 | 4.71e-3 | 1.15e-8 | 1.06e-9 | 7.38e-8 |
 | sqrt | 1.77e-7 | 1.16e-7 | 4.74e-7 | 2.70e-12 | 1.78e-12 | 7.16e-12 |
 <!-- ACCURACY_END -->

src/bounded.rs

@@ -77,7 +77,7 @@ impl<T: CordicNumber> NonNegative<T> {
         Self(x_sq.saturating_sub(T::one()))
     }
 
-    /// Returns the inner value.
+    /// Unwraps the inner value.
     #[inline]
     #[must_use]
     pub const fn get(self) -> T {
@@ -103,7 +103,7 @@ impl<T: CordicNumber> UnitInterval<T> {
         (value >= -one && value <= one).then_some(Self(value))
     }
 
-    /// Returns the inner value.
+    /// Unwraps the inner value.
     #[inline]
     #[must_use]
     pub const fn get(self) -> T {
@@ -158,13 +158,13 @@ impl<T: CordicNumber> OpenUnitInterval<T> {
     /// `(x - 1) / (x + 1)` is in `(-1/3, 1/3)` which is a subset of `(-1, 1)`.
     #[inline]
     #[must_use]
-    pub fn from_normalized_ln_arg(x: NormalizedLnArg<T>) -> Self {
+    pub(crate) fn from_normalized_ln_arg(x: NormalizedLnArg<T>) -> Self {
         let x_minus_1 = x.0 - T::one();
         let x_plus_1 = x.0 + T::one();
         Self(x_minus_1.div(x_plus_1))
     }
 
-    /// Returns the inner value.
+    /// Unwraps the inner value.
     #[inline]
     #[must_use]
     pub const fn get(self) -> T {
@@ -188,7 +188,7 @@ impl<T: CordicNumber> AtLeastOne<T> {
         (value >= T::one()).then_some(Self(value))
     }
 
-    /// Returns the inner value.
+    /// Unwraps the inner value.
     #[inline]
     #[must_use]
     pub const fn get(self) -> T {
@@ -201,25 +201,18 @@ impl<T: CordicNumber> AtLeastOne<T> {
 /// After normalizing the input for ln computation, the value is always
 /// in this range, which guarantees that `(x-1)/(x+1)` is in `(-1/3, 1/3)`.
 #[derive(Clone, Copy, Debug)]
-pub struct NormalizedLnArg<T>(T);
+pub(crate) struct NormalizedLnArg<T>(T);
 
 impl<T: CordicNumber> NormalizedLnArg<T> {
     /// Creates a new `NormalizedLnArg` from the normalization loop result.
     ///
     /// The ln function's normalization loop guarantees the result is in [0.5, 2].
-    /// This constructor trusts that invariant (used only in ln implementation).
+    /// This constructor trusts that invariant.
     #[inline]
     #[must_use]
     pub(crate) const fn from_normalized(value: T) -> Self {
         Self(value)
     }
-
-    /// Returns the inner value.
-    #[inline]
-    #[must_use]
-    pub const fn get(self) -> T {
-        self.0
-    }
 }
 
 #[cfg(test)]
@@ -326,12 +319,6 @@ mod tests {
         assert_eq!(at_least.get(), I16F16::from_num(2));
     }
 
-    #[test]
-    fn normalized_ln_arg_get() {
-        let norm = NormalizedLnArg::from_normalized(I16F16::from_num(1.5));
-        assert_eq!(norm.get(), I16F16::from_num(1.5));
-    }
-
     #[test]
     fn open_unit_interval_from_normalized_ln_arg() {
         let norm = NormalizedLnArg::from_normalized(I16F16::from_num(1.5));

src/error.rs

@@ -15,12 +15,6 @@ pub enum Error {
         /// Human-readable description of the valid domain.
         expected: &'static str,
     },
-
-    /// The computation would overflow the fixed-point representation.
-    Overflow {
-        /// Name of the function that encountered the error.
-        function: &'static str,
-    },
 }
 
 impl Error {
@@ -29,12 +23,6 @@ impl Error {
     pub const fn domain(function: &'static str, expected: &'static str) -> Self {
         Self::DomainError { function, expected }
     }
-
-    /// Create an overflow error for the given function.
-    #[must_use]
-    pub const fn overflow(function: &'static str) -> Self {
-        Self::Overflow { function }
-    }
 }
 
 impl fmt::Display for Error {
@@ -46,9 +34,6 @@ impl fmt::Display for Error {
                     "{function}: input outside valid domain, expected {expected}"
                 )
             }
-            Self::Overflow { function } => {
-                write!(f, "{function}: result would overflow")
-            }
         }
     }
 }

src/kernel/cordic.rs

@@ -1,12 +1,11 @@
 //! Core CORDIC iteration implementations.
 //!
-//! The CORDIC algorithm operates in three modes, each with two directions:
+//! The CORDIC algorithm operates in two modes, each with two directions:
 //!
 //! | Mode | Rotation (z → 0) | Vectoring (y → 0) |
 //! |------|------------------|-------------------|
 //! | Circular | sin, cos | atan |
 //! | Hyperbolic | sinh, cosh | atanh, ln |
-//! | Linear | multiply | divide |
 //!
 //! # Algorithm
 //!
@@ -30,16 +29,13 @@ use crate::traits::CordicNumber;
 
 /// Table lookup for CORDIC iteration.
 ///
-/// # Safety Invariant
 /// 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.
-/// The bounds check is optimized away in release builds.
 #[inline]
 const fn table_lookup(table: &[i64; 64], index: u32) -> i64 {
-    // SAFETY: Iteration limits guarantee index < 64 (see doc comment above).
     #[allow(
         clippy::indexing_slicing,
         reason = "index bounded by CORDIC iteration limits"
@@ -109,13 +105,11 @@ pub fn circular_rotation<T: CordicNumber>(mut x: T, mut y: T, mut z: T) -> (T, T
         let angle = T::from_i1f63(table_lookup(&ATAN_TABLE, i));
 
         if z >= zero {
-            // Rotate counter-clockwise (positive angle)
             let x_new = x.saturating_sub(y >> i);
             y = y.saturating_add(x >> i);
             x = x_new;
             z -= angle;
         } else {
-            // Rotate clockwise (negative angle)
             let x_new = x.saturating_add(y >> i);
             y = y.saturating_sub(x >> i);
             x = x_new;

src/ops/circular.rs

@@ -12,20 +12,28 @@ pub fn sin_cos<T: CordicNumber>(angle: T) -> (T, T) {
     let pi = T::pi();
     let frac_pi_2 = T::frac_pi_2();
     let zero = T::zero();
-
-    // Reduce angle to [-π, π] range first.
-    let mut reduced = angle;
     let two_pi = pi + pi;
-    let mut i = 0;
-    while reduced > pi && i < 64 {
-        reduced -= two_pi;
-        i += 1;
-    }
-    i = 0;
-    while reduced < -pi && i < 64 {
-        reduced += two_pi;
-        i += 1;
-    }
+
+    // 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);
+        let n = quotient.round();
+        angle.saturating_sub(n.saturating_mul(two_pi))
+    } else {
+        angle
+    };
+
+    // Clamp to [-π, π] to handle any residual from saturation.
+    // This is a safety net; mathematically unnecessary for valid inputs.
+    let reduced = if reduced > pi {
+        reduced.saturating_sub(two_pi)
+    } else if reduced < -pi {
+        reduced.saturating_add(two_pi)
+    } else {
+        reduced
+    };
 
     // Further reduce to [-π/2, π/2] and track sign
     let (reduced, negate) = if reduced > frac_pi_2 {
@@ -61,7 +69,44 @@ pub fn cos<T: CordicNumber>(angle: T) -> T {
     sin_cos(angle).1
 }
 
-/// Tangent. May overflow near ±π/2.
+/// Tangent. Returns `sin(angle) / cos(angle)`.
+///
+/// # Overflow Behavior
+///
+/// Tangent has poles at ±π/2, ±3π/2, etc. where it approaches ±∞.
+/// Since these poles occur at irrational values that cannot be exactly
+/// represented in fixed-point, this function will never divide by
+/// exactly zero. However, near poles the result may:
+///
+/// - Saturate to `T::MAX` or `T::MIN` for very small denominators
+/// - Produce very large finite values that may overflow in subsequent operations
+///
+/// The threshold for potential overflow is approximately:
+/// - `|angle - π/2| < 2^(-frac_bits/2)` for the nearest pole
+///
+/// For I16F16, this means angles within ~0.004 radians of π/2 may overflow.
+/// For I32F32, within ~0.00003 radians.
+///
+/// If you need to detect near-pole conditions, check `cos(angle).abs()`
+/// against a threshold before calling `tan`.
+///
+/// # Example
+///
+/// ```
+/// use fixed::types::I16F16;
+/// use fixed_analytics::{tan, cos};
+///
+/// let angle = I16F16::from_num(1.5); // Close to π/2 ≈ 1.571
+///
+/// // Safe: check cosine magnitude first
+/// let c = cos(angle);
+/// if c.abs() > I16F16::from_num(0.01) {
+///     let t = tan(angle);
+///     // Use t safely
+/// } else {
+///     // Handle near-pole case
+/// }
+/// ```
 #[must_use]
 pub fn tan<T: CordicNumber>(angle: T) -> T {
     let (s, c) = sin_cos(angle);