From 118f44a3141d3c725ab667262fa6b109e16a8c09 Mon Sep 17 00:00:00 2001
From: Joshua Simmons <josh@nega.tv>
Date: Sat, 15 Oct 2022 23:34:23 +0200
Subject: [PATCH] Fix some edge case handling in `sin_cos_pi`

Update tests as well.
---
 narcissus-maths/src/sin_cos_pi.rs       | 16 +++++---
 narcissus-maths/tests/exhaustive_f32.rs | 54 +++++++++++++++++--------
 2 files changed, 49 insertions(+), 21 deletions(-)

diff --git a/narcissus-maths/src/sin_cos_pi.rs b/narcissus-maths/src/sin_cos_pi.rs
index 9d6f5fd..ce6de9f 100644
--- a/narcissus-maths/src/sin_cos_pi.rs
+++ b/narcissus-maths/src/sin_cos_pi.rs
@@ -42,6 +42,9 @@ pub fn sin_cos_pi_f32(a: f32) -> (f32, f32) {
     const S: [f32; 4] = F32_SIN_PI_7_K;
     const C: [f32; 4] = F32_COS_PI_8_K;
 
+    // cos_pi(a) = 1.0f for |a| > 2^24, but cos_pi(Inf) = NaN
+    let a = if a.abs() < 16777216.0 { a } else { a * 0.0 };
+
     // Range reduction.
     let r = (a + a).round();
     let i: u32 = unsafe { r.to_int_unchecked() };
@@ -67,11 +70,12 @@ pub fn sin_cos_pi_f32(a: f32) -> (f32, f32) {
     let s = s.mul_add(r2, S[1]);
     let s = r.mul_add(S[0], r3 * s);
 
-    if i & 1 != 0 {
-        (c, s)
-    } else {
-        (s, c)
-    }
+    let (s, c) = if i & 1 != 0 { (c, s) } else { (s, c) };
+
+    // IEEE-754: sin_pi(+n) is +0 and sin_pi(-n) is -0 for positive integers n
+    let s = if a == a.floor() { a * 0.0 } else { s };
+
+    (s, c)
 }
 
 #[cfg(test)]
@@ -80,6 +84,8 @@ mod tests {
 
     #[test]
     fn basics() {
+        assert_eq!(sin_cos_pi_f32(f32::from_bits(0x7f4135c6)), (0.0, 1.0));
+
         assert_eq!(sin_cos_pi_f32(-1.5), (1.0, 0.0));
         assert_eq!(sin_cos_pi_f32(-1.0), (0.0, -1.0));
         assert_eq!(sin_cos_pi_f32(-0.5), (-1.0, 0.0));
diff --git a/narcissus-maths/tests/exhaustive_f32.rs b/narcissus-maths/tests/exhaustive_f32.rs
index 9d9f85c..02295f6 100644
--- a/narcissus-maths/tests/exhaustive_f32.rs
+++ b/narcissus-maths/tests/exhaustive_f32.rs
@@ -77,6 +77,7 @@ struct FloatErrors {
 }
 
 const PREC: u32 = 280;
+const COUNT: u32 = 2_139_095_040_u32;
 
 struct FloatCheck {
     reference: fn(&mut rug::Float),
@@ -124,8 +125,6 @@ impl FloatCheck {
         let ref_val = self.tmp.to_f32_round(rug::float::Round::Nearest);
 
         if our_val != ref_val && !(our_val.is_nan() && ref_val.is_nan()) {
-            println!("u: {u:#08x}, our_val: {our_val}, ref_val: {ref_val}");
-
             self.errors.num_errors += 1;
             let e = ulp_error_f32(our_val, ref_val, x, self.reference);
             if e > self.errors.max_error_ulp {
@@ -136,9 +135,12 @@ impl FloatCheck {
     }
 }
 
-fn check_exhaustive_f32(reference: fn(&mut rug::Float), ours: fn(f32) -> f32, tan_mode: bool) {
-    const COUNT: u32 = 2_139_095_040_u32;
-    const SPLIT: u32 = 256;
+fn check_exhaustive_f32(
+    reference: fn(&mut rug::Float),
+    ours: fn(f32) -> f32,
+    tan_mode: bool,
+) -> FloatErrors {
+    const SPLIT: u32 = 512;
 
     // Generate ranges.
     assert_eq!(COUNT % SPLIT, 0);
@@ -205,14 +207,7 @@ fn check_exhaustive_f32(reference: fn(&mut rug::Float), ours: fn(f32) -> f32, ta
         }
     });
 
-    println!(
-        "errors: {} ({:.2}%)",
-        errors.num_errors,
-        (errors.num_errors as f64 / (COUNT * 2) as f64) * 100.0
-    );
-    println!("max error ulps: {}", errors.max_error_ulp);
-    println!("max error error at: {:#08x}", errors.max_error_val);
-    assert_eq!(errors.num_errors, 0);
+    errors
 }
 
 fn ref_sin_pi_f32(x: &mut rug::Float) {
@@ -230,17 +225,44 @@ fn ref_tan_pi_f32(x: &mut rug::Float) {
 #[test]
 #[ignore]
 pub fn exhaustive_sin_pi() {
-    check_exhaustive_f32(ref_sin_pi_f32, |a| sin_cos_pi_f32(a).0, false)
+    let errors = check_exhaustive_f32(ref_sin_pi_f32, |a| sin_cos_pi_f32(a).0, false);
+    println!(
+        "errors: {} ({:.2}%)",
+        errors.num_errors,
+        (errors.num_errors as f64 / (COUNT * 2) as f64) * 100.0
+    );
+    println!("max error ulps: {}", errors.max_error_ulp);
+    println!("max error error at: {:#08x}", errors.max_error_val);
+    assert_eq!(errors.max_error_ulp, 1);
+    assert_eq!(errors.num_errors, 715_786_430);
 }
 
 #[test]
 #[ignore]
 pub fn exhaustive_cos_pi() {
-    check_exhaustive_f32(ref_cos_pi_f32, |a| sin_cos_pi_f32(a).1, false)
+    let errors = check_exhaustive_f32(ref_cos_pi_f32, |a| sin_cos_pi_f32(a).1, false);
+    println!(
+        "errors: {} ({:.2}%)",
+        errors.num_errors,
+        (errors.num_errors as f64 / (COUNT * 2) as f64) * 100.0
+    );
+    println!("max error ulps: {}", errors.max_error_ulp);
+    println!("max error error at: {:#08x}", errors.max_error_val);
+    assert_eq!(errors.num_errors, 33_455_772);
+    assert_eq!(errors.max_error_ulp, 1);
 }
 
 #[test]
 #[ignore]
 pub fn exhaustive_tan_pi() {
-    check_exhaustive_f32(ref_tan_pi_f32, tan_pi_f32, true)
+    let errors = check_exhaustive_f32(ref_tan_pi_f32, tan_pi_f32, true);
+    println!(
+        "errors: {} ({:.2}%)",
+        errors.num_errors,
+        (errors.num_errors as f64 / (COUNT * 2) as f64) * 100.0
+    );
+    println!("max error ulps: {}", errors.max_error_ulp);
+    println!("max error error at: {:#08x}", errors.max_error_val);
+    assert_eq!(errors.num_errors, 68_630_324);
+    assert_eq!(errors.max_error_ulp, 2);
 }
-- 
2.49.0