Skip to content

v2.0.0: Replace CORDIC rotation with polynomial evaluation #7

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
GeEom merged 2 commits into from
Mar 16, 2026
Merged

v2.0.0: Replace CORDIC rotation with polynomial evaluation #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c47a8b6
Replace CORDIC rotation with polynomial evaluation
GeEom Mar 16, 2026
5ff2e53
Linting
GeEom Mar 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Cargo.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[package]
name = "fixed_analytics"
version = "1.0.0"
version = "2.0.0"
edition = "2024"
rust-version = "1.88"
authors = ["David Gathercole"]
Expand Down
32 changes: 16 additions & 16 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,14 +30,14 @@ Requires Rust 1.88 or later.

```toml
[dependencies]
fixed_analytics = "1.0.0"
fixed_analytics = "2.0.0"
```

For `no_std` environments:

```toml
[dependencies]
fixed_analytics = { version = "1.0.0", default-features = false }
fixed_analytics = { version = "2.0.0", default-features = false }
```

## Available Functions
Expand All @@ -54,21 +54,21 @@ fixed_analytics = { version = "1.0.0", default-features = false }
| Exponential | `exp`, `pow2` | `ln`, `log2`, `log10` |
| Algebraic | — | `sqrt` |

Functions are calculated via CORDIC, Newton-Raphson, and Taylor series techniques. Complete absence of panic is verified at the linker level via the [`no-panic`](https://github.com/dtolnay/no-panic) crate.
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.

### Saturation Behavior

The following total functions saturate, clamping to the representable range near the following thresholds.

| Function | I16F16 Threshold | I32F32 Threshold | Result |
|----------|------------------|------------------|--------|
| `exp` | x ≥ 22.2 | x ≥ 44.4 | `T::MAX` |
| `exp` | x ≤ -9.2 | x ≤ -16.2 | Zero |
| `exp` | x ≥ 10.4 | x ≥ 21.5 | `T::MAX` |
| `exp` | x ≤ -11.1 | x ≤ -22.2 | Zero |
| `pow2` | x ≥ 15.0 | x ≥ 31.0 | `T::MAX` |
| `pow2` | x ≤ -13.2 | x ≤ -23.3 | Zero |
| `pow2` | x ≤ -16.1 | x ≤ -32.1 | 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` |
| `tan` | \|x - pole\| < 4e-5 | \|x - pole\| < 5e-10 | `T::MAX` or `T::MIN` |

Where for `tan`, "pole" refers to ±π/2, ±3π/2, ±5π/2, ...

Expand All @@ -79,24 +79,24 @@ Relative error statistics measured against MPFR reference implementations. Accur

| Function | I16F16 Mean | I16F16 Median | I16F16 P95 | I32F32 Mean | I32F32 Median | I32F32 P95 |
|----------|-------------|---------------|------------|-------------|---------------|------------|
| sin | 6.19e-4 | 9.34e-5 | 1.30e-3 | 1.22e-8 | 1.79e-9 | 2.52e-8 |
| cos | 6.82e-4 | 1.02e-4 | 1.46e-3 | 1.30e-8 | 1.91e-9 | 2.83e-8 |
| tan | 2.47e-4 | 9.84e-5 | 6.70e-4 | 6.00e-9 | 1.89e-9 | 1.40e-8 |
| 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 | 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.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 |
| 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 |
| exp | 1.14e-2 | 6.67e-5 | 7.87e-2 | 1.91e-7 | 2.24e-9 | 1.30e-6 |
| exp | 1.14e-2 | 2.32e-5 | 7.88e-2 | 1.91e-7 | 1.73e-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.28e-4 | 4.74e-5 | 4.71e-3 | 1.15e-8 | 1.06e-9 | 7.38e-8 |
| pow2 | 7.21e-4 | 2.58e-5 | 4.72e-3 | 1.13e-8 | 4.93e-10 | 7.43e-8 |
| sqrt | 1.77e-7 | 1.16e-7 | 4.74e-7 | 2.70e-12 | 1.78e-12 | 7.16e-12 |
<!-- ACCURACY_END -->
161 changes: 6 additions & 155 deletions src/kernel/cordic.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
//! Core CORDIC iteration implementations.
//!
//! The CORDIC algorithm operates in two modes, each with two directions:
//! CORDIC vectoring mode drives y toward zero while accumulating angles:
//!
//! | Mode | Rotation (z → 0) | Vectoring (y → 0) |
//! |------|------------------|-------------------|
//! | Circular | sin, cos | atan |
//! | Hyperbolic | sinh, cosh | atanh, ln |
//! | Mode | Vectoring (y → 0) |
//! |------|-------------------|
//! | Circular | atan |
//! | Hyperbolic | atanh, ln |
//!
//! # Algorithm
//!
Expand All @@ -22,9 +22,7 @@
//! - angle[i] = atan(2^-i) for circular, atanh(2^-i) for hyperbolic

use crate::tables::hyperbolic::needs_repeat;
use crate::tables::{
ATAN_TABLE, ATANH_TABLE, CIRCULAR_GAIN_INV, HYPERBOLIC_GAIN, HYPERBOLIC_GAIN_INV,
};
use crate::tables::{ATAN_TABLE, ATANH_TABLE};
use crate::traits::CordicNumber;

/// Table lookup for CORDIC iteration.
Expand All @@ -43,83 +41,6 @@ const fn table_lookup(table: &[i64; 64], index: u32) -> i64 {
table[index as usize]
}

/// Returns the CORDIC scale factor (1/K ≈ 0.6073).
///
/// Pre-multiply initial vectors by this to compensate for CORDIC gain.
#[inline]
#[must_use]
pub fn cordic_scale_factor<T: CordicNumber>() -> T {
T::from_i1f63(CIRCULAR_GAIN_INV)
}

/// Returns the hyperbolic gain factor (`K_h` ≈ 0.8282).
///
/// After hyperbolic CORDIC iterations, results are scaled by `1/K_h`.
/// To compensate, divide by `K_h` (or multiply by `1/K_h`).
#[inline]
#[must_use]
pub fn hyperbolic_gain<T: CordicNumber>() -> T {
T::from_i1f63(HYPERBOLIC_GAIN)
}

/// Returns the inverse hyperbolic gain factor (`1/K_h` ≈ 1.2075).
///
/// Pre-multiply initial vectors by this to compensate for hyperbolic CORDIC gain.
/// This uses a precomputed constant, avoiding runtime division.
#[inline]
#[must_use]
pub fn hyperbolic_gain_inv<T: CordicNumber>() -> T {
T::from_i2f62(HYPERBOLIC_GAIN_INV)
}

/// Performs circular CORDIC in rotation mode.
///
/// Given an initial vector (x, y) and angle z, rotates the vector by angle z.
/// After iteration:
/// - x ≈ K * (x₀ * cos(z₀) - y₀ * sin(z₀))
/// - y ≈ K * (y₀ * cos(z₀) + x₀ * sin(z₀))
/// - z ≈ 0
///
/// Where K is the circular gain factor (~1.6468).
///
/// # Arguments
///
/// * `x` - Initial x coordinate
/// * `y` - Initial y coordinate
/// * `z` - Angle to rotate by (in radians)
///
/// # Returns
///
/// Tuple of (x, y, z) after CORDIC iterations.
///
/// # Note
///
/// The input angle should be in the range [-1.74, 1.74] radians for
/// convergence. Use argument reduction for larger angles.
#[must_use]
pub fn circular_rotation<T: CordicNumber>(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 {
let angle = T::from_i1f63(table_lookup(&ATAN_TABLE, i));

if z >= zero {
let x_new = x.saturating_sub(y >> i);
y = y.saturating_add(x >> i);
x = x_new;
z -= angle;
} else {
let x_new = x.saturating_add(y >> i);
y = y.saturating_sub(x >> i);
x = x_new;
z += angle;
}
}

(x, y, z)
}

/// Performs circular CORDIC in vectoring mode.
///
/// Given an initial vector (x, y), rotates it until y ≈ 0.
Expand Down Expand Up @@ -167,76 +88,6 @@ pub fn circular_vectoring<T: CordicNumber>(mut x: T, mut y: T, mut z: T) -> (T,
(x, y, z)
}

/// Performs hyperbolic CORDIC in rotation mode.
///
/// Given initial values (x, y, z), performs hyperbolic pseudo-rotations
/// to drive z toward zero.
///
/// After iteration:
/// - x ≈ `K_h` * (x₀ * cosh(z₀) + y₀ * sinh(z₀))
/// - y ≈ `K_h` * (y₀ * cosh(z₀) + x₀ * sinh(z₀))
/// - z ≈ 0
///
/// Where `K_h` is the hyperbolic gain factor (~1.2075).
///
/// # Arguments
///
/// * `x` - Initial x value
/// * `y` - Initial y value
/// * `z` - Hyperbolic angle to "rotate" by
///
/// # Returns
///
/// Tuple of (x, y, z) after CORDIC iterations.
///
/// # Note
///
/// - Hyperbolic CORDIC starts at i=1 (not i=0)
/// - Certain iterations must be repeated for convergence
/// - Input z should be in range [-1.12, 1.12] for convergence
#[must_use]
pub fn hyperbolic_rotation<T: CordicNumber>(mut x: T, mut y: T, mut z: T) -> (T, T, T) {
let zero = T::zero();
// Use frac_bits iterations, capped at 54 for table bounds.
let max_iterations = T::frac_bits().min(54);

let mut i: u32 = 1; // Hyperbolic starts at i=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));

if z >= zero {
// "Rotate" in positive direction
let x_new = x.saturating_add(y >> i);
y = y.saturating_add(x >> i);
x = x_new;
z -= angle;
} else {
// "Rotate" in negative direction
let x_new = x.saturating_sub(y >> i);
y = y.saturating_sub(x >> i);
x = x_new;
z += angle;
}

iteration_count += 1;

// Handle repetition for convergence
if needs_repeat(i) && !repeated {
repeated = true;
// Don't increment i, repeat this iteration
} else {
repeated = false;
i += 1;
}
}

(x, y, z)
}

/// Performs hyperbolic CORDIC in vectoring mode.
///
/// Drives y toward zero while accumulating the hyperbolic angle.
Expand Down
16 changes: 6 additions & 10 deletions src/kernel/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,22 +10,18 @@
//! z' = z - σ·atan(2^-i)
//! ```
//!
//! **Rotation mode** (z→0): computes sin/cos from angle.
//! **Vectoring mode** (y→0): computes atan from coordinates.
//! **Vectoring mode** (y→0): computes atan/atanh from coordinates.
//!
//! Hyperbolic mode uses `atanh(2^-i)` tables and requires iteration repeats
//! at indices 4, 13, 40, ... for convergence.
//!
//! | Mode | Rotation (z → 0) | Vectoring (y → 0) |
//! |------|------------------|-------------------|
//! | Circular | sin, cos | atan |
//! | Hyperbolic | sinh, cosh | atanh, ln |
//! | Mode | Vectoring (y → 0) |
//! |------|-------------------|
//! | Circular | atan |
//! | Hyperbolic | atanh, ln |
//!
//! Users should call functions in [`crate::ops`] rather than kernels directly.

mod cordic;

pub use crate::kernel::cordic::{circular_rotation, circular_vectoring, cordic_scale_factor};
pub use crate::kernel::cordic::{
hyperbolic_gain, hyperbolic_gain_inv, hyperbolic_rotation, hyperbolic_vectoring,
};
pub use crate::kernel::cordic::{circular_vectoring, hyperbolic_vectoring};
53 changes: 48 additions & 5 deletions src/ops/circular.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,9 @@

use crate::bounded::{NonNegative, UnitInterval};
use crate::error::{Error, Result};
use crate::kernel::{circular_rotation, circular_vectoring, cordic_scale_factor};
use crate::kernel::circular_vectoring;
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;

/// Sine and cosine. More efficient than separate calls. Accepts any angle.
Expand All @@ -12,7 +13,6 @@ use crate::traits::CordicNumber;
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();
let two_pi = pi + pi;

// Reduce angle to [-π, π] using direct quotient computation.
Expand Down Expand Up @@ -45,9 +45,52 @@ pub fn sin_cos<T: CordicNumber>(angle: T) -> (T, T) {
(reduced, false)
};

// Run CORDIC with unit vector scaled by inverse gain
let inv_gain = cordic_scale_factor();
let (cos_val, sin_val, _) = circular_rotation(inv_gain, zero, reduced);
// Polynomial evaluation via factored Horner form.
// To avoid catastrophic cancellation near π/2, reduce to [0, π/4]:
// For |x| ∈ [0, π/4]: sin(x) = sin_poly(x), cos(x) = cos_poly(x)
// For |x| ∈ (π/4, π/2]: sin(x) = cos_poly(π/2-|x|), cos(x) = sin_poly(π/2-|x|)
let one = T::one();
let frac_pi_4 = T::frac_pi_4();
let abs_reduced = reduced.abs();
let (poly_arg, swapped) = if abs_reduced >= frac_pi_4 {
(frac_pi_2.saturating_sub(abs_reduced), true)
} else {
(abs_reduced, false)
};
let u = poly_arg.saturating_mul(poly_arg);

// Evaluate sin and cos polynomials over [0, π/4] using minimax
// (Chebyshev) coefficients. Uses multiply-by-constant instead of
// division, avoiding cumulative rounding error from per-step divides.
//
// sin(x) = x + x³·P(x²) where P = minimax poly of (sin(x)-x)/x³
// cos(x) = 1 + x²·Q(x²) where Q = minimax poly of (cos(x)-1)/x²
let (sp_val, cp_val) = if T::frac_bits() >= 24 {
// High precision: degree 15 sin, degree 14 cos
let sp = horner(&SIN_P_HI, u);
let sin_approx = poly_arg.saturating_add(poly_arg.saturating_mul(u).saturating_mul(sp));
let cp = horner(&COS_Q_HI, u);
(sin_approx, one.saturating_add(u.saturating_mul(cp)))
} else {
// Low precision: degree 9 sin, degree 8 cos
let sp = horner(&SIN_P_LO, u);
let sin_approx = poly_arg.saturating_add(poly_arg.saturating_mul(u).saturating_mul(sp));
let cp = horner(&COS_Q_LO, u);
(sin_approx, one.saturating_add(u.saturating_mul(cp)))
};

// Map back: if we swapped, sin(x) = cos_poly, cos(x) = sin_poly
// Also restore sign of sin for negative angles.
let (sin_unsigned, cos_val) = if swapped {
(cp_val, sp_val)
} else {
(sp_val, cp_val)
};
let sin_val = if reduced < T::zero() {
-sin_unsigned
} else {
sin_unsigned
};

if negate {
(-sin_val, -cos_val)
Expand Down
Loading
Loading