Skip to content

Commit 6952d0f

Browse files
authored
Merge pull request #7 from GeEom/further_accuracy
v2.0.0: Replace CORDIC rotation with polynomial evaluation
2 parents 4d13c1b + 5ff2e53 commit 6952d0f

22 files changed

Lines changed: 586 additions & 649 deletions

Cargo.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
[package]
22
name = "fixed_analytics"
3-
version = "1.0.0"
3+
version = "2.0.0"
44
edition = "2024"
55
rust-version = "1.88"
66
authors = ["David Gathercole"]

README.md

Lines changed: 16 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -30,14 +30,14 @@ Requires Rust 1.88 or later.
3030

3131
```toml
3232
[dependencies]
33-
fixed_analytics = "1.0.0"
33+
fixed_analytics = "2.0.0"
3434
```
3535

3636
For `no_std` environments:
3737

3838
```toml
3939
[dependencies]
40-
fixed_analytics = { version = "1.0.0", default-features = false }
40+
fixed_analytics = { version = "2.0.0", default-features = false }
4141
```
4242

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

57-
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.
57+
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.
5858

5959
### Saturation Behavior
6060

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

6363
| Function | I16F16 Threshold | I32F32 Threshold | Result |
6464
|----------|------------------|------------------|--------|
65-
| `exp` | x ≥ 22.2 | x ≥ 44.4 | `T::MAX` |
66-
| `exp` | x ≤ -9.2 | x ≤ -16.2 | Zero |
65+
| `exp` | x ≥ 10.4 | x ≥ 21.5 | `T::MAX` |
66+
| `exp` | x ≤ -11.1 | x ≤ -22.2 | Zero |
6767
| `pow2` | x ≥ 15.0 | x ≥ 31.0 | `T::MAX` |
68-
| `pow2` | x ≤ -13.2 | x ≤ -23.3 | Zero |
68+
| `pow2` | x ≤ -16.1 | x ≤ -32.1 | Zero |
6969
| `sinh` | \|x\| ≥ 11.1 | \|x\| ≥ 22.2 | `T::MAX` or `T::MIN` |
7070
| `cosh` | \|x\| ≥ 11.1 | \|x\| ≥ 22.2 | `T::MAX` |
71-
| `tan` | \|x - pole\| < 8e-5 | \|x - pole\| < 5e-10 | `T::MAX` or `T::MIN` |
71+
| `tan` | \|x - pole\| < 4e-5 | \|x - pole\| < 5e-10 | `T::MAX` or `T::MIN` |
7272

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

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

8080
| Function | I16F16 Mean | I16F16 Median | I16F16 P95 | I32F32 Mean | I32F32 Median | I32F32 P95 |
8181
|----------|-------------|---------------|------------|-------------|---------------|------------|
82-
| sin | 6.19e-4 | 9.34e-5 | 1.30e-3 | 1.22e-8 | 1.79e-9 | 2.52e-8 |
83-
| cos | 6.82e-4 | 1.02e-4 | 1.46e-3 | 1.30e-8 | 1.91e-9 | 2.83e-8 |
84-
| tan | 2.47e-4 | 9.84e-5 | 6.70e-4 | 6.00e-9 | 1.89e-9 | 1.40e-8 |
82+
| sin | 6.06e-4 | 8.78e-5 | 1.28e-3 | 1.16e-8 | 1.68e-9 | 2.43e-8 |
83+
| cos | 6.45e-4 | 9.03e-5 | 1.38e-3 | 1.22e-8 | 1.72e-9 | 2.64e-8 |
84+
| tan | 7.20e-5 | 3.57e-5 | 2.20e-4 | 1.28e-9 | 3.98e-10 | 3.03e-9 |
8585
| asin | 2.87e-4 | 5.93e-5 | 6.46e-4 | 5.34e-9 | 8.82e-10 | 1.03e-8 |
8686
| acos | 3.61e-5 | 2.18e-5 | 1.14e-4 | 5.37e-10 | 3.19e-10 | 1.71e-9 |
8787
| atan | 2.71e-5 | 2.21e-5 | 6.29e-5 | 3.69e-10 | 2.92e-10 | 8.74e-10 |
88-
| sinh | 1.82e-4 | 1.34e-4 | 5.03e-4 | 1.05e-8 | 2.35e-9 | 9.30e-9 |
89-
| cosh | 1.73e-4 | 1.23e-4 | 5.00e-4 | 1.02e-8 | 2.07e-9 | 9.16e-9 |
90-
| tanh | 2.08e-5 | 1.38e-5 | 5.89e-5 | 1.64e-9 | 1.31e-10 | 1.23e-9 |
91-
| coth | 1.20e-5 | 4.83e-6 | 5.57e-5 | 4.00e-10 | 1.39e-10 | 1.30e-9 |
88+
| sinh | 9.80e-5 | 6.23e-5 | 2.79e-4 | 1.52e-9 | 9.64e-10 | 4.29e-9 |
89+
| cosh | 9.40e-5 | 5.75e-5 | 2.77e-4 | 1.44e-9 | 8.90e-10 | 4.25e-9 |
90+
| tanh | 1.60e-5 | 1.32e-5 | 2.56e-5 | 2.25e-10 | 1.22e-10 | 3.90e-10 |
91+
| coth | 6.68e-6 | 3.54e-6 | 1.80e-5 | 1.41e-10 | 1.16e-10 | 2.74e-10 |
9292
| asinh | 6.44e-4 | 4.83e-4 | 1.75e-3 | 1.03e-8 | 7.59e-9 | 2.85e-8 |
9393
| acosh | 6.74e-4 | 5.21e-4 | 1.80e-3 | 1.05e-8 | 7.96e-9 | 2.88e-8 |
9494
| atanh | 3.01e-4 | 5.90e-5 | 6.25e-4 | 6.68e-9 | 1.32e-9 | 1.44e-8 |
9595
| acoth | 2.10e-3 | 1.33e-3 | 6.67e-3 | 4.26e-8 | 2.62e-8 | 1.39e-7 |
96-
| exp | 1.14e-2 | 6.67e-5 | 7.87e-2 | 1.91e-7 | 2.24e-9 | 1.30e-6 |
96+
| exp | 1.14e-2 | 2.32e-5 | 7.88e-2 | 1.91e-7 | 1.73e-9 | 1.30e-6 |
9797
| ln | 1.35e-5 | 8.76e-6 | 2.97e-5 | 4.50e-10 | 3.48e-10 | 9.17e-10 |
9898
| log2 | 1.33e-5 | 8.48e-6 | 2.92e-5 | 3.46e-10 | 2.24e-10 | 7.21e-10 |
9999
| log10 | 1.44e-5 | 9.28e-6 | 3.14e-5 | 4.49e-10 | 3.27e-10 | 9.07e-10 |
100-
| pow2 | 7.28e-4 | 4.74e-5 | 4.71e-3 | 1.15e-8 | 1.06e-9 | 7.38e-8 |
100+
| pow2 | 7.21e-4 | 2.58e-5 | 4.72e-3 | 1.13e-8 | 4.93e-10 | 7.43e-8 |
101101
| sqrt | 1.77e-7 | 1.16e-7 | 4.74e-7 | 2.70e-12 | 1.78e-12 | 7.16e-12 |
102102
<!-- ACCURACY_END -->

src/kernel/cordic.rs

Lines changed: 6 additions & 155 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,11 @@
11
//! Core CORDIC iteration implementations.
22
//!
3-
//! The CORDIC algorithm operates in two modes, each with two directions:
3+
//! CORDIC vectoring mode drives y toward zero while accumulating angles:
44
//!
5-
//! | Mode | Rotation (z → 0) | Vectoring (y → 0) |
6-
//! |------|------------------|-------------------|
7-
//! | Circular | sin, cos | atan |
8-
//! | Hyperbolic | sinh, cosh | atanh, ln |
5+
//! | Mode | Vectoring (y → 0) |
6+
//! |------|-------------------|
7+
//! | Circular | atan |
8+
//! | Hyperbolic | atanh, ln |
99
//!
1010
//! # Algorithm
1111
//!
@@ -22,9 +22,7 @@
2222
//! - angle[i] = atan(2^-i) for circular, atanh(2^-i) for hyperbolic
2323
2424
use crate::tables::hyperbolic::needs_repeat;
25-
use crate::tables::{
26-
ATAN_TABLE, ATANH_TABLE, CIRCULAR_GAIN_INV, HYPERBOLIC_GAIN, HYPERBOLIC_GAIN_INV,
27-
};
25+
use crate::tables::{ATAN_TABLE, ATANH_TABLE};
2826
use crate::traits::CordicNumber;
2927

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

46-
/// Returns the CORDIC scale factor (1/K ≈ 0.6073).
47-
///
48-
/// Pre-multiply initial vectors by this to compensate for CORDIC gain.
49-
#[inline]
50-
#[must_use]
51-
pub fn cordic_scale_factor<T: CordicNumber>() -> T {
52-
T::from_i1f63(CIRCULAR_GAIN_INV)
53-
}
54-
55-
/// Returns the hyperbolic gain factor (`K_h` ≈ 0.8282).
56-
///
57-
/// After hyperbolic CORDIC iterations, results are scaled by `1/K_h`.
58-
/// To compensate, divide by `K_h` (or multiply by `1/K_h`).
59-
#[inline]
60-
#[must_use]
61-
pub fn hyperbolic_gain<T: CordicNumber>() -> T {
62-
T::from_i1f63(HYPERBOLIC_GAIN)
63-
}
64-
65-
/// Returns the inverse hyperbolic gain factor (`1/K_h` ≈ 1.2075).
66-
///
67-
/// Pre-multiply initial vectors by this to compensate for hyperbolic CORDIC gain.
68-
/// This uses a precomputed constant, avoiding runtime division.
69-
#[inline]
70-
#[must_use]
71-
pub fn hyperbolic_gain_inv<T: CordicNumber>() -> T {
72-
T::from_i2f62(HYPERBOLIC_GAIN_INV)
73-
}
74-
75-
/// Performs circular CORDIC in rotation mode.
76-
///
77-
/// Given an initial vector (x, y) and angle z, rotates the vector by angle z.
78-
/// After iteration:
79-
/// - x ≈ K * (x₀ * cos(z₀) - y₀ * sin(z₀))
80-
/// - y ≈ K * (y₀ * cos(z₀) + x₀ * sin(z₀))
81-
/// - z ≈ 0
82-
///
83-
/// Where K is the circular gain factor (~1.6468).
84-
///
85-
/// # Arguments
86-
///
87-
/// * `x` - Initial x coordinate
88-
/// * `y` - Initial y coordinate
89-
/// * `z` - Angle to rotate by (in radians)
90-
///
91-
/// # Returns
92-
///
93-
/// Tuple of (x, y, z) after CORDIC iterations.
94-
///
95-
/// # Note
96-
///
97-
/// The input angle should be in the range [-1.74, 1.74] radians for
98-
/// convergence. Use argument reduction for larger angles.
99-
#[must_use]
100-
pub fn circular_rotation<T: CordicNumber>(mut x: T, mut y: T, mut z: T) -> (T, T, T) {
101-
let zero = T::zero();
102-
let iterations = T::frac_bits().min(62);
103-
104-
for i in 0..iterations {
105-
let angle = T::from_i1f63(table_lookup(&ATAN_TABLE, i));
106-
107-
if z >= zero {
108-
let x_new = x.saturating_sub(y >> i);
109-
y = y.saturating_add(x >> i);
110-
x = x_new;
111-
z -= angle;
112-
} else {
113-
let x_new = x.saturating_add(y >> i);
114-
y = y.saturating_sub(x >> i);
115-
x = x_new;
116-
z += angle;
117-
}
118-
}
119-
120-
(x, y, z)
121-
}
122-
12344
/// Performs circular CORDIC in vectoring mode.
12445
///
12546
/// Given an initial vector (x, y), rotates it until y ≈ 0.
@@ -167,76 +88,6 @@ pub fn circular_vectoring<T: CordicNumber>(mut x: T, mut y: T, mut z: T) -> (T,
16788
(x, y, z)
16889
}
16990

170-
/// Performs hyperbolic CORDIC in rotation mode.
171-
///
172-
/// Given initial values (x, y, z), performs hyperbolic pseudo-rotations
173-
/// to drive z toward zero.
174-
///
175-
/// After iteration:
176-
/// - x ≈ `K_h` * (x₀ * cosh(z₀) + y₀ * sinh(z₀))
177-
/// - y ≈ `K_h` * (y₀ * cosh(z₀) + x₀ * sinh(z₀))
178-
/// - z ≈ 0
179-
///
180-
/// Where `K_h` is the hyperbolic gain factor (~1.2075).
181-
///
182-
/// # Arguments
183-
///
184-
/// * `x` - Initial x value
185-
/// * `y` - Initial y value
186-
/// * `z` - Hyperbolic angle to "rotate" by
187-
///
188-
/// # Returns
189-
///
190-
/// Tuple of (x, y, z) after CORDIC iterations.
191-
///
192-
/// # Note
193-
///
194-
/// - Hyperbolic CORDIC starts at i=1 (not i=0)
195-
/// - Certain iterations must be repeated for convergence
196-
/// - Input z should be in range [-1.12, 1.12] for convergence
197-
#[must_use]
198-
pub fn hyperbolic_rotation<T: CordicNumber>(mut x: T, mut y: T, mut z: T) -> (T, T, T) {
199-
let zero = T::zero();
200-
// Use frac_bits iterations, capped at 54 for table bounds.
201-
let max_iterations = T::frac_bits().min(54);
202-
203-
let mut i: u32 = 1; // Hyperbolic starts at i=1
204-
let mut iteration_count: u32 = 0;
205-
let mut repeated = false;
206-
207-
while iteration_count < max_iterations && i < 64 {
208-
let table_index = i.saturating_sub(1);
209-
let angle = T::from_i1f63(table_lookup(&ATANH_TABLE, table_index));
210-
211-
if z >= zero {
212-
// "Rotate" in positive direction
213-
let x_new = x.saturating_add(y >> i);
214-
y = y.saturating_add(x >> i);
215-
x = x_new;
216-
z -= angle;
217-
} else {
218-
// "Rotate" in negative direction
219-
let x_new = x.saturating_sub(y >> i);
220-
y = y.saturating_sub(x >> i);
221-
x = x_new;
222-
z += angle;
223-
}
224-
225-
iteration_count += 1;
226-
227-
// Handle repetition for convergence
228-
if needs_repeat(i) && !repeated {
229-
repeated = true;
230-
// Don't increment i, repeat this iteration
231-
} else {
232-
repeated = false;
233-
i += 1;
234-
}
235-
}
236-
237-
(x, y, z)
238-
}
239-
24091
/// Performs hyperbolic CORDIC in vectoring mode.
24192
///
24293
/// Drives y toward zero while accumulating the hyperbolic angle.

src/kernel/mod.rs

Lines changed: 6 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -10,22 +10,18 @@
1010
//! z' = z - σ·atan(2^-i)
1111
//! ```
1212
//!
13-
//! **Rotation mode** (z→0): computes sin/cos from angle.
14-
//! **Vectoring mode** (y→0): computes atan from coordinates.
13+
//! **Vectoring mode** (y→0): computes atan/atanh from coordinates.
1514
//!
1615
//! Hyperbolic mode uses `atanh(2^-i)` tables and requires iteration repeats
1716
//! at indices 4, 13, 40, ... for convergence.
1817
//!
19-
//! | Mode | Rotation (z → 0) | Vectoring (y → 0) |
20-
//! |------|------------------|-------------------|
21-
//! | Circular | sin, cos | atan |
22-
//! | Hyperbolic | sinh, cosh | atanh, ln |
18+
//! | Mode | Vectoring (y → 0) |
19+
//! |------|-------------------|
20+
//! | Circular | atan |
21+
//! | Hyperbolic | atanh, ln |
2322
//!
2423
//! Users should call functions in [`crate::ops`] rather than kernels directly.
2524
2625
mod cordic;
2726

28-
pub use crate::kernel::cordic::{circular_rotation, circular_vectoring, cordic_scale_factor};
29-
pub use crate::kernel::cordic::{
30-
hyperbolic_gain, hyperbolic_gain_inv, hyperbolic_rotation, hyperbolic_vectoring,
31-
};
27+
pub use crate::kernel::cordic::{circular_vectoring, hyperbolic_vectoring};

src/ops/circular.rs

Lines changed: 48 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,9 @@
22
33
use crate::bounded::{NonNegative, UnitInterval};
44
use crate::error::{Error, Result};
5-
use crate::kernel::{circular_rotation, circular_vectoring, cordic_scale_factor};
5+
use crate::kernel::circular_vectoring;
66
use crate::ops::algebraic::sqrt_nonneg;
7+
use crate::tables::chebyshev::{COS_Q_HI, COS_Q_LO, SIN_P_HI, SIN_P_LO, horner};
78
use crate::traits::CordicNumber;
89

910
/// Sine and cosine. More efficient than separate calls. Accepts any angle.
@@ -12,7 +13,6 @@ use crate::traits::CordicNumber;
1213
pub fn sin_cos<T: CordicNumber>(angle: T) -> (T, T) {
1314
let pi = T::pi();
1415
let frac_pi_2 = T::frac_pi_2();
15-
let zero = T::zero();
1616
let two_pi = pi + pi;
1717

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

48-
// Run CORDIC with unit vector scaled by inverse gain
49-
let inv_gain = cordic_scale_factor();
50-
let (cos_val, sin_val, _) = circular_rotation(inv_gain, zero, reduced);
48+
// Polynomial evaluation via factored Horner form.
49+
// To avoid catastrophic cancellation near π/2, reduce to [0, π/4]:
50+
// For |x| ∈ [0, π/4]: sin(x) = sin_poly(x), cos(x) = cos_poly(x)
51+
// For |x| ∈ (π/4, π/2]: sin(x) = cos_poly(π/2-|x|), cos(x) = sin_poly(π/2-|x|)
52+
let one = T::one();
53+
let frac_pi_4 = T::frac_pi_4();
54+
let abs_reduced = reduced.abs();
55+
let (poly_arg, swapped) = if abs_reduced >= frac_pi_4 {
56+
(frac_pi_2.saturating_sub(abs_reduced), true)
57+
} else {
58+
(abs_reduced, false)
59+
};
60+
let u = poly_arg.saturating_mul(poly_arg);
61+
62+
// Evaluate sin and cos polynomials over [0, π/4] using minimax
63+
// (Chebyshev) coefficients. Uses multiply-by-constant instead of
64+
// division, avoiding cumulative rounding error from per-step divides.
65+
//
66+
// sin(x) = x + x³·P(x²) where P = minimax poly of (sin(x)-x)/x³
67+
// cos(x) = 1 + x²·Q(x²) where Q = minimax poly of (cos(x)-1)/x²
68+
let (sp_val, cp_val) = if T::frac_bits() >= 24 {
69+
// High precision: degree 15 sin, degree 14 cos
70+
let sp = horner(&SIN_P_HI, u);
71+
let sin_approx = poly_arg.saturating_add(poly_arg.saturating_mul(u).saturating_mul(sp));
72+
let cp = horner(&COS_Q_HI, u);
73+
(sin_approx, one.saturating_add(u.saturating_mul(cp)))
74+
} else {
75+
// Low precision: degree 9 sin, degree 8 cos
76+
let sp = horner(&SIN_P_LO, u);
77+
let sin_approx = poly_arg.saturating_add(poly_arg.saturating_mul(u).saturating_mul(sp));
78+
let cp = horner(&COS_Q_LO, u);
79+
(sin_approx, one.saturating_add(u.saturating_mul(cp)))
80+
};
81+
82+
// Map back: if we swapped, sin(x) = cos_poly, cos(x) = sin_poly
83+
// Also restore sign of sin for negative angles.
84+
let (sin_unsigned, cos_val) = if swapped {
85+
(cp_val, sp_val)
86+
} else {
87+
(sp_val, cp_val)
88+
};
89+
let sin_val = if reduced < T::zero() {
90+
-sin_unsigned
91+
} else {
92+
sin_unsigned
93+
};
5194

5295
if negate {
5396
(-sin_val, -cos_val)

0 commit comments

Comments
 (0)