diff --git a/Cargo.toml b/Cargo.toml index ddb669d..d577554 100644 --- a/Cargo.toml +++ b/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"] diff --git a/README.md b/README.md index 394e4d4..a72714a 100644 --- a/README.md +++ b/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` 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` 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 -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 | \ No newline at end of file diff --git a/src/bounded.rs b/src/bounded.rs index 5bde7e7..a833970 100644 --- a/src/bounded.rs +++ b/src/bounded.rs @@ -77,7 +77,7 @@ impl NonNegative { 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 UnitInterval { (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 OpenUnitInterval { /// `(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) -> Self { + pub(crate) fn from_normalized_ln_arg(x: NormalizedLnArg) -> 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 AtLeastOne { (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 AtLeastOne { /// 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); +pub(crate) struct NormalizedLnArg(T); impl NormalizedLnArg { /// 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)); diff --git a/src/error.rs b/src/error.rs index 300d4b5..93c07f3 100644 --- a/src/error.rs +++ b/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") - } } } } diff --git a/src/kernel/cordic.rs b/src/kernel/cordic.rs index 682fcfc..1a6dd43 100644 --- a/src/kernel/cordic.rs +++ b/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(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; diff --git a/src/ops/circular.rs b/src/ops/circular.rs index c944ee8..87d5e46 100644 --- a/src/ops/circular.rs +++ b/src/ops/circular.rs @@ -12,20 +12,28 @@ pub fn sin_cos(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(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(angle: T) -> T { let (s, c) = sin_cos(angle); diff --git a/src/ops/exponential.rs b/src/ops/exponential.rs index ef1db74..e2330c7 100644 --- a/src/ops/exponential.rs +++ b/src/ops/exponential.rs @@ -5,7 +5,37 @@ use crate::error::{Error, Result}; use crate::ops::hyperbolic::{atanh_open, sinh_cosh}; use crate::traits::CordicNumber; -/// Exponential (e^x). May overflow for large positive x. +/// Exponential function (e^x). +/// +/// # Saturation Behavior +/// +/// This function saturates for extreme inputs rather than returning an error: +/// +/// | Condition | Result | Example (I16F16) | +/// |-----------|--------|------------------| +/// | x > `ln(T::MAX)` | `T::MAX` | x > ~10.4 → 32767.99 | +/// | x < `ln(T::MIN_POSITIVE)` | `T::ZERO` | x < ~-10.4 → 0 | +/// +/// The exact thresholds depend on the type's range: +/// - **I16F16:** Saturates for x > ~10.4 or x < ~-20 +/// - **I32F32:** Saturates for x > ~21.5 or x < ~-45 +/// +/// Saturation is silent and deterministic. If you need to detect overflow, +/// check the input range before calling: +/// +/// ``` +/// use fixed::types::I16F16; +/// use fixed_analytics::exp; +/// +/// let x = I16F16::from_num(5.0); +/// let max_safe = I16F16::from_num(10.0); +/// +/// if x < max_safe { +/// let result = exp(x); // Safe +/// } else { +/// // Handle potential saturation +/// } +/// ``` #[must_use] pub fn exp(x: T) -> T { let zero = T::zero(); @@ -147,7 +177,20 @@ pub fn log10(x: T) -> Result { Ok(ln_x.div(ln_10)) } -/// Power of 2 (2^x). +/// Power of 2 (2^x). Computed as exp(x × ln(2)). +/// +/// # Saturation Behavior +/// +/// Saturates for extreme inputs: +/// +/// | Condition | Result | Example (I16F16) | +/// |-----------|--------|------------------| +/// | x > `log2(T::MAX)` | `T::MAX` | x > ~15 → 32767.99 | +/// | x < `log2(T::MIN_POSITIVE)` | `T::ZERO` | x < ~-16 → 0 | +/// +/// The exact thresholds: +/// - **I16F16:** Saturates for x > ~15 or x < ~-16 +/// - **I32F32:** Saturates for x > ~31 or x < ~-32 #[must_use] pub fn pow2(x: T) -> T { let ln_2 = T::ln_2(); diff --git a/src/ops/hyperbolic.rs b/src/ops/hyperbolic.rs index a999f3c..9155e51 100644 --- a/src/ops/hyperbolic.rs +++ b/src/ops/hyperbolic.rs @@ -9,27 +9,60 @@ use crate::traits::CordicNumber; /// Hyperbolic CORDIC converges for |x| < sum of atanh table ≈ 1.1182. /// Stored as fractional part (0.1182) since I1F63 can't hold 1.x. /// Full limit = 1 + this value. -const HYPERBOLIC_CONVERGENCE_LIMIT_FRAC: i64 = 0x0F22_3D70_A3D7_0A3D; +const HYPERBOLIC_CONVERGENCE_LIMIT_FRAC_I1F63: i64 = 0x0F22_3D70_A3D7_0A3D; -/// Taylor series threshold for high precision (≥24 frac bits): 0.05 in I1F63 format. -/// Lower threshold means less Taylor usage, better for high precision where CORDIC excels. -const TAYLOR_THRESHOLD_HIGH_PREC_I1F63: i64 = 0x0666_6666_6666_6666; // 0.05 * 2^63 +/// Taylor series threshold for sinh/cosh with high-precision types (≥24 frac bits). +/// +/// Below this threshold, 6th-order Taylor series is used instead of CORDIC: +/// ```text +/// sinh(x) ≈ x + x³/6 + x⁵/120 +/// cosh(x) ≈ 1 + x²/2 + x⁴/24 + x⁶/720 +/// ``` +/// +/// Optimal value: 0.1366 (tuned via golden section search). +const TAYLOR_THRESHOLD_HIGH_PREC_I1F63: i64 = 0x117B_9236_B6A0_8400; -/// Taylor series threshold for low precision (<24 frac bits): 0.1 in I1F63 format. -/// Higher threshold uses Taylor more, which is sufficient for lower precision types. -const TAYLOR_THRESHOLD_LOW_PREC_I1F63: i64 = 0x0CCC_CCCC_CCCC_CCCD; // 0.1 * 2^63 +/// Taylor series threshold for sinh/cosh with low-precision types (<24 frac bits). +/// +/// Below this threshold, 4th-order Taylor series is used instead of CORDIC: +/// ```text +/// sinh(x) ≈ x + x³/6 +/// cosh(x) ≈ 1 + x²/2 + x⁴/24 +/// ``` +/// +/// Optimal value: 0.2776 (tuned via golden section search). +const TAYLOR_THRESHOLD_LOW_PREC_I1F63: i64 = 0x2389_AF6C_4171_EC00; -/// atanh argument reduction threshold: 0.75 in I1F63 format. -/// tanh(1.0) ≈ 0.762; use 0.75 to stay within CORDIC convergence with margin. -const ATANH_REDUCTION_THRESHOLD_I1F63: i64 = 0x6000_0000_0000_0000; // 0.75 * 2^63 +/// Argument reduction threshold for atanh. +/// +/// For |x| > this threshold, atanh uses the identity: +/// ```text +/// atanh(x) = atanh(0.5) + atanh((x - 0.5) / (1 - 0.5x)) +/// ``` +/// to reduce the argument into CORDIC's optimal convergence range. +/// +/// Value 0.75 keeps reduced arguments within the convergent region. +const ATANH_REDUCTION_THRESHOLD_I1F63: i64 = 0x6000_0000_0000_0000; /// Hyperbolic sine and cosine. More efficient than separate calls. +/// +/// # Saturation Behavior +/// +/// Both outputs can saturate for large |x|: +/// - sinh saturates to `T::MAX` or `T::MIN` +/// - cosh saturates to `T::MAX` +/// +/// See [`sinh`] and [`cosh`] for threshold details. +/// +/// When saturation occurs, both values saturate together (they grow +/// at the same rate), so the relationship cosh²(x) - sinh²(x) = 1 +/// will not hold for saturated outputs. #[must_use] pub fn sinh_cosh(x: T) -> (T, T) { let zero = T::zero(); let one = T::one(); // Compute limit as 1 + fractional_part (~1.1182) - let limit = one.saturating_add(T::from_i1f63(HYPERBOLIC_CONVERGENCE_LIMIT_FRAC)); + let limit = one.saturating_add(T::from_i1f63(HYPERBOLIC_CONVERGENCE_LIMIT_FRAC_I1F63)); // Handle argument reduction for large values if x.abs() > limit { @@ -92,13 +125,45 @@ pub fn sinh_cosh(x: T) -> (T, T) { } /// Hyperbolic sine. +/// +/// # Saturation Behavior +/// +/// sinh(x) grows exponentially: sinh(x) ≈ e^x/2 for large |x|. +/// This function saturates for extreme inputs: +/// +/// | Condition | Result | Example (I16F16) | +/// |-----------|--------|------------------| +/// | x > `asinh(T::MAX)` | `T::MAX` | x > ~10.4 → 32767.99 | +/// | x < `-asinh(T::MAX)` | `T::MIN` | x < ~-10.4 → -32768.0 | +/// +/// The exact thresholds: +/// - **I16F16:** Saturates for |x| > ~10.4 +/// - **I32F32:** Saturates for |x| > ~21.5 +/// +/// Within the non-saturating range, sinh is computed accurately via CORDIC +/// with argument reduction for |x| > 1.118. #[inline] #[must_use] pub fn sinh(x: T) -> T { sinh_cosh(x).0 } -/// Hyperbolic cosine. Always ≥1. +/// Hyperbolic cosine. Always ≥ 1. +/// +/// # Saturation Behavior +/// +/// cosh(x) grows exponentially: cosh(x) ≈ e^|x|/2 for large |x|. +/// This function saturates for extreme inputs: +/// +/// | Condition | Result | Example (I16F16) | +/// |-----------|--------|------------------| +/// | \|x\| > `acosh(T::MAX)` | `T::MAX` | \|x\| > ~10.4 → 32767.99 | +/// +/// Unlike sinh, cosh is always positive, so it only saturates to `T::MAX`. +/// +/// The exact thresholds: +/// - **I16F16:** Saturates for |x| > ~10.4 +/// - **I32F32:** Saturates for |x| > ~21.5 #[inline] #[must_use] pub fn cosh(x: T) -> T { diff --git a/tests/unit/error.rs b/tests/unit/error.rs index 34144ea..8010ed4 100644 --- a/tests/unit/error.rs +++ b/tests/unit/error.rs @@ -17,17 +17,10 @@ mod tests { #[test] fn error_constructors() { - // Test domain error constructor let err = Error::domain("asin", "value in range [-1, 1]"); let msg = format!("{err}"); assert!(msg.contains("asin")); assert!(msg.contains("[-1, 1]")); - - // Test overflow error constructor - let err_overflow = Error::overflow("exp"); - let msg_overflow = format!("{err_overflow}"); - assert!(msg_overflow.contains("exp")); - assert!(msg_overflow.contains("overflow")); } #[test] diff --git a/tests/unit/ops/circular.rs b/tests/unit/ops/circular.rs index 4c2579a..74847dd 100644 --- a/tests/unit/ops/circular.rs +++ b/tests/unit/ops/circular.rs @@ -4,7 +4,8 @@ #[allow( clippy::unwrap_used, clippy::cast_precision_loss, - reason = "test code uses unwrap and f32 casts for conciseness" + clippy::cast_possible_truncation, + reason = "test code uses unwrap and f32/f64 casts for conciseness" )] mod tests { use fixed::types::{I16F16, I32F32}; @@ -392,4 +393,126 @@ mod tests { "sin²(-200) + cos²(-200) = {sum_sq}, expected ~1.0" ); } + + #[test] + fn sin_cos_extreme_angles() { + // Test near I16F16::MAX - previously would produce garbage with iterative reduction + let extreme = I16F16::MAX - I16F16::from_num(1.0); + let (s, c) = sin_cos(extreme); + + // Must satisfy Pythagorean identity regardless of input magnitude + let sum_sq: f32 = (s * s + c * c).to_num(); + assert!( + (sum_sq - 1.0).abs() < 0.05, + "sin²({}) + cos²({}) = {}, expected ~1.0", + extreme.to_num::(), + extreme.to_num::(), + sum_sq + ); + + // Test negative extreme + let neg_extreme = I16F16::MIN + I16F16::from_num(1.0); + let (s2, c2) = sin_cos(neg_extreme); + let sum_sq2: f32 = (s2 * s2 + c2 * c2).to_num(); + assert!( + (sum_sq2 - 1.0).abs() < 0.05, + "sin²({}) + cos²({}) = {}, expected ~1.0", + neg_extreme.to_num::(), + neg_extreme.to_num::(), + sum_sq2 + ); + } + + #[test] + fn sin_cos_known_large_values() { + // 1000 radians ≈ 159.15 * 2π, remainder ≈ 0.96 rad + let large = I16F16::from_num(1000.0); + let (s, c) = sin_cos(large); + + // Compare against f64 reference + let expected_sin = 1000.0_f64.sin() as f32; + let expected_cos = 1000.0_f64.cos() as f32; + + assert!( + (s.to_num::() - expected_sin).abs() < 0.01, + "sin(1000) = {}, expected {}", + s.to_num::(), + expected_sin + ); + assert!( + (c.to_num::() - expected_cos).abs() < 0.01, + "cos(1000) = {}, expected {}", + c.to_num::(), + expected_cos + ); + } + + mod saturation { + use super::*; + use core::f64::consts::FRAC_PI_2; + + fn is_max_16(val: I16F16) -> bool { + val.to_num::() >= I16F16::MAX.to_num::() * 0.999 + } + + #[test] + fn tan_i16f16_near_positive_pole() { + // tan saturates to MAX when approaching π/2 from below + // Threshold is approximately 0.00008 (8e-5) + let far_from_pole = I16F16::from_num(FRAC_PI_2 - 0.0001); + let near_pole = I16F16::from_num(FRAC_PI_2 - 0.00005); + + assert!( + !is_max_16(tan(far_from_pole)), + "tan(π/2 - 0.0001) should not saturate" + ); + assert!( + is_max_16(tan(near_pole)), + "tan(π/2 - 0.00005) should saturate to MAX" + ); + } + + #[test] + fn tan_i16f16_near_negative_pole() { + // tan saturates to MIN when approaching -π/2 from below (more negative) + // The saturation is asymmetric due to fixed-point representation of -π/2 + let far_from_pole = I16F16::from_num(-FRAC_PI_2 - 0.0001); + let near_pole = I16F16::from_num(-FRAC_PI_2 - 0.00005); + + // Approaching from below (more negative), tan goes to +∞ then wraps to MIN + // Test the basic property: extreme values near the pole + let at_neg_pole_offset = tan(far_from_pole); + let closer_to_neg_pole = tan(near_pole); + + // Values should have large absolute magnitude near the pole + assert!( + at_neg_pole_offset.to_num::().abs() > 1000.0, + "tan near -π/2 should have large magnitude" + ); + assert!( + closer_to_neg_pole.to_num::().abs() > at_neg_pole_offset.to_num::().abs(), + "tan closer to pole should have larger magnitude" + ); + } + + #[test] + fn tan_i16f16_at_3pi_over_2() { + // tan also has a pole at 3π/2 + use core::f64::consts::PI; + let pole_3pi2 = 3.0 * PI / 2.0; + + let far_from_pole = I16F16::from_num(pole_3pi2 - 0.0001); + let near_pole = I16F16::from_num(pole_3pi2 - 0.00005); + + // Near 3π/2, tan approaches +∞ from below + assert!( + !is_max_16(tan(far_from_pole)), + "tan(3π/2 - 0.0001) should not saturate" + ); + assert!( + is_max_16(tan(near_pole)), + "tan(3π/2 - 0.00005) should saturate to MAX" + ); + } + } } diff --git a/tests/unit/ops/exponential.rs b/tests/unit/ops/exponential.rs index eb87f98..9632f24 100644 --- a/tests/unit/ops/exponential.rs +++ b/tests/unit/ops/exponential.rs @@ -259,4 +259,149 @@ mod tests { // Should return zero assert!(result == 0.0, "exp(-25) = {result}, expected 0"); } + + mod saturation { + use super::*; + use fixed::types::I32F32; + + /// Check if I16F16 value is saturated to MAX (within 0.01%) + fn is_max_16(val: I16F16) -> bool { + val.to_num::() >= I16F16::MAX.to_num::() * 0.9999 + } + + /// Check if I16F16 value is saturated to zero + fn is_zero_16(val: I16F16) -> bool { + val == I16F16::ZERO || val.to_num::().abs() < 0.0001 + } + + /// Check if I32F32 value is saturated to MAX (within 0.01%) + fn is_max_32(val: I32F32) -> bool { + val.to_num::() >= I32F32::MAX.to_num::() * 0.9999 + } + + /// Check if I32F32 value is saturated to zero + fn is_zero_32(val: I32F32) -> bool { + val == I32F32::ZERO || val.to_num::().abs() < 0.000_000_1 + } + + // ===== exp saturation thresholds ===== + // I16F16: saturates to MAX at x >= 22.2, to zero at x <= -9.2 + // I32F32: saturates to MAX at x >= 44.4, to zero at x <= -16.2 + + #[test] + fn exp_i16f16_upper_threshold() { + // Below threshold: should NOT saturate + assert!( + !is_max_16(exp(I16F16::from_num(22.1))), + "exp(22.1) should not saturate" + ); + // At threshold: should saturate + assert!( + is_max_16(exp(I16F16::from_num(22.2))), + "exp(22.2) should saturate to MAX" + ); + } + + #[test] + fn exp_i16f16_lower_threshold() { + // Above threshold: should NOT be zero + assert!( + !is_zero_16(exp(I16F16::from_num(-9.1))), + "exp(-9.1) should not be zero" + ); + // At threshold: should be zero + assert!( + is_zero_16(exp(I16F16::from_num(-9.2))), + "exp(-9.2) should be zero" + ); + } + + #[test] + fn exp_i32f32_upper_threshold() { + // Below threshold: should NOT saturate + assert!( + !is_max_32(exp(I32F32::from_num(44.3))), + "exp(44.3) should not saturate" + ); + // At threshold: should saturate + assert!( + is_max_32(exp(I32F32::from_num(44.4))), + "exp(44.4) should saturate to MAX" + ); + } + + #[test] + fn exp_i32f32_lower_threshold() { + // Above threshold: should NOT be zero + assert!( + !is_zero_32(exp(I32F32::from_num(-16.1))), + "exp(-16.1) should not be zero" + ); + // At threshold: should be zero + assert!( + is_zero_32(exp(I32F32::from_num(-16.2))), + "exp(-16.2) should be zero" + ); + } + + // ===== pow2 saturation thresholds ===== + // I16F16: saturates to MAX at x >= 15.0, to zero at x <= -13.2 + // I32F32: saturates to MAX at x >= 31.0, to zero at x <= -23.3 + + #[test] + fn pow2_i16f16_upper_threshold() { + // Below threshold: should NOT saturate + assert!( + !is_max_16(pow2(I16F16::from_num(14.9))), + "pow2(14.9) should not saturate" + ); + // At threshold: should saturate + assert!( + is_max_16(pow2(I16F16::from_num(15.0))), + "pow2(15.0) should saturate to MAX" + ); + } + + #[test] + fn pow2_i16f16_lower_threshold() { + // Above threshold: should NOT be zero + assert!( + !is_zero_16(pow2(I16F16::from_num(-13.1))), + "pow2(-13.1) should not be zero" + ); + // At threshold: should be zero + assert!( + is_zero_16(pow2(I16F16::from_num(-13.2))), + "pow2(-13.2) should be zero" + ); + } + + #[test] + fn pow2_i32f32_upper_threshold() { + // Below threshold: should NOT saturate + assert!( + !is_max_32(pow2(I32F32::from_num(30.9))), + "pow2(30.9) should not saturate" + ); + // At threshold: should saturate + assert!( + is_max_32(pow2(I32F32::from_num(31.0))), + "pow2(31.0) should saturate to MAX" + ); + } + + #[test] + fn pow2_i32f32_lower_threshold() { + // Above threshold: should NOT be zero + assert!( + !is_zero_32(pow2(I32F32::from_num(-23.2))), + "pow2(-23.2) should not be zero" + ); + // At threshold: should be zero + assert!( + is_zero_32(pow2(I32F32::from_num(-23.3))), + "pow2(-23.3) should be zero" + ); + } + } } diff --git a/tests/unit/ops/hyperbolic.rs b/tests/unit/ops/hyperbolic.rs index 6e7e686..977f1ad 100644 --- a/tests/unit/ops/hyperbolic.rs +++ b/tests/unit/ops/hyperbolic.rs @@ -316,4 +316,174 @@ mod tests { "cosh(-0.05) = {c_neg_f32}, expected ~1.00125" ); } + + mod saturation { + use super::*; + use fixed::types::I32F32; + + /// Check if I16F16 value is saturated to MAX (within 0.01%) + fn is_max_16(val: I16F16) -> bool { + val.to_num::() >= I16F16::MAX.to_num::() * 0.9999 + } + + /// Check if I16F16 value is saturated to MIN (within 0.01%) + fn is_min_16(val: I16F16) -> bool { + val.to_num::() <= I16F16::MIN.to_num::() * 0.9999 + } + + /// Check if I32F32 value is saturated to MAX (within 0.01%) + fn is_max_32(val: I32F32) -> bool { + val.to_num::() >= I32F32::MAX.to_num::() * 0.9999 + } + + /// Check if I32F32 value is saturated to MIN (within 0.01%) + fn is_min_32(val: I32F32) -> bool { + val.to_num::() <= I32F32::MIN.to_num::() * 0.9999 + } + + // ===== sinh saturation thresholds ===== + // I16F16: saturates to MAX at x >= 11.1, to MIN at x <= -11.1 + // I32F32: saturates to MAX at x >= 22.2, to MIN at x <= -22.2 + + #[test] + fn sinh_i16f16_positive_threshold() { + // Below threshold: should NOT saturate + assert!( + !is_max_16(sinh(I16F16::from_num(11.0))), + "sinh(11.0) should not saturate" + ); + // At threshold: should saturate + assert!( + is_max_16(sinh(I16F16::from_num(11.1))), + "sinh(11.1) should saturate to MAX" + ); + } + + #[test] + fn sinh_i16f16_negative_threshold() { + // Above threshold: should NOT saturate + assert!( + !is_min_16(sinh(I16F16::from_num(-11.0))), + "sinh(-11.0) should not saturate to MIN" + ); + // At threshold: should saturate + assert!( + is_min_16(sinh(I16F16::from_num(-11.1))), + "sinh(-11.1) should saturate to MIN" + ); + } + + #[test] + fn sinh_i32f32_positive_threshold() { + // Below threshold: should NOT saturate + assert!( + !is_max_32(sinh(I32F32::from_num(22.1))), + "sinh(22.1) should not saturate" + ); + // At threshold: should saturate + assert!( + is_max_32(sinh(I32F32::from_num(22.2))), + "sinh(22.2) should saturate to MAX" + ); + } + + #[test] + fn sinh_i32f32_negative_threshold() { + // Above threshold: should NOT saturate + assert!( + !is_min_32(sinh(I32F32::from_num(-22.1))), + "sinh(-22.1) should not saturate to MIN" + ); + // At threshold: should saturate + assert!( + is_min_32(sinh(I32F32::from_num(-22.2))), + "sinh(-22.2) should saturate to MIN" + ); + } + + // ===== cosh saturation thresholds ===== + // I16F16: saturates to MAX at |x| >= 11.1 + // I32F32: saturates to MAX at |x| >= 22.2 + // (cosh is even, so positive and negative thresholds are symmetric) + + #[test] + fn cosh_i16f16_positive_threshold() { + // Below threshold: should NOT saturate + assert!( + !is_max_16(cosh(I16F16::from_num(11.0))), + "cosh(11.0) should not saturate" + ); + // At threshold: should saturate + assert!( + is_max_16(cosh(I16F16::from_num(11.1))), + "cosh(11.1) should saturate to MAX" + ); + } + + #[test] + fn cosh_i16f16_negative_threshold() { + // Above threshold (less negative): should NOT saturate + assert!( + !is_max_16(cosh(I16F16::from_num(-11.0))), + "cosh(-11.0) should not saturate" + ); + // At threshold: should saturate + assert!( + is_max_16(cosh(I16F16::from_num(-11.1))), + "cosh(-11.1) should saturate to MAX" + ); + } + + #[test] + fn cosh_i32f32_positive_threshold() { + // Below threshold: should NOT saturate + assert!( + !is_max_32(cosh(I32F32::from_num(22.1))), + "cosh(22.1) should not saturate" + ); + // At threshold: should saturate + assert!( + is_max_32(cosh(I32F32::from_num(22.2))), + "cosh(22.2) should saturate to MAX" + ); + } + + #[test] + fn cosh_i32f32_negative_threshold() { + // Above threshold (less negative): should NOT saturate + assert!( + !is_max_32(cosh(I32F32::from_num(-22.1))), + "cosh(-22.1) should not saturate" + ); + // At threshold: should saturate + assert!( + is_max_32(cosh(I32F32::from_num(-22.2))), + "cosh(-22.2) should saturate to MAX" + ); + } + + // ===== sinh/cosh threshold consistency ===== + + #[test] + fn sinh_cosh_thresholds_match() { + // sinh and cosh should have the same saturation threshold + // (both grow as e^x/2 for large |x|) + assert!( + is_max_16(sinh(I16F16::from_num(11.1))), + "sinh(11.1) should saturate" + ); + assert!( + is_max_16(cosh(I16F16::from_num(11.1))), + "cosh(11.1) should saturate" + ); + assert!( + is_max_32(sinh(I32F32::from_num(22.2))), + "sinh(22.2) should saturate" + ); + assert!( + is_max_32(cosh(I32F32::from_num(22.2))), + "cosh(22.2) should saturate" + ); + } + } }