From: Joshua Simmons <josh@nega.tv>
Date: Sat, 22 Oct 2022 14:25:44 +0000 (+0200)
Subject: Improve error rate on sin_pi
X-Git-Url: https://git.nega.tv//gitweb.cgi?a=commitdiff_plain;h=fe09be0d551a34875ac8ce7f5a814cbd9dc795bb;p=josh%2Fnarcissus

Improve error rate on sin_pi

Similarly to our tan implementation, use a 48 bit leading coefficient.
Reduces the percentage of faithfully rounded results from 17.4% to 1.3%.
---

diff --git a/narcissus-maths/doc/sincostan.sollya b/narcissus-maths/doc/sincostan.sollya
index eb3d517..ef8280b 100644
--- a/narcissus-maths/doc/sincostan.sollya
+++ b/narcissus-maths/doc/sincostan.sollya
@@ -121,6 +121,8 @@ procedure f32_data_write(name,a,f,r)
 
 // Actual start of building polynomial
 
+B = [|D,SG...|];
+
 print("sin(pi x) in binary32");
 I = [|1,3,5,7|];
 F = sin(pi*x);
@@ -129,6 +131,8 @@ f32_data_write("F32_SIN_PI_7", P,F,R);
 f32_write_list("F32_SIN_PI_7", list_of_odd_coeff(P));
 print("\n");
 
+B = [|SG...|];
+
 print("cos(pi x) in binary32");
 I = [|0,2,4,6,8|];
 F = cos(pi*x);
diff --git a/narcissus-maths/src/sin_cos_pi.rs b/narcissus-maths/src/sin_cos_pi.rs
index 0f1450f..e9912cb 100644
--- a/narcissus-maths/src/sin_cos_pi.rs
+++ b/narcissus-maths/src/sin_cos_pi.rs
@@ -3,12 +3,11 @@
 // Sollya code for generating these polynomials is in `doc/sincostan.sollya`
 
 // constants for sin(pi x), cos(pi x) for x on [-1/4,1/4]
-const F32_SIN_PI_7_K: [f32; 4] = unsafe {
-    std::mem::transmute::<[u32; 4], _>([
-        0x40490fdb, // 0x1.921fb6p1
-        0xc0a55df6, // -0x1.4abbecp2
-        0x40233590, // 0x1.466b2p1
-        0xbf17acc9, // -0x1.2f5992p-1
+const F32_SIN_PI_7_K: [f32; 3] = unsafe {
+    std::mem::transmute::<[u32; 3], _>([
+        0xc0a55ddd, // -0x1.4abbbap2
+        0x40232f6e, // 0x1.465edcp1
+        0xbf16cf2d, // -0x1.2d9e5ap-1
     ])
 };
 
@@ -42,7 +41,7 @@ fn mulsign_f32(x: f32, s: u32) -> f32 {
 /// assert_eq!(cos, 1.0);
 /// ```
 pub fn sin_cos_pi_f32(a: f32) -> (f32, f32) {
-    const S: [f32; 4] = F32_SIN_PI_7_K;
+    const S: [f32; 3] = 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
@@ -67,11 +66,10 @@ pub fn sin_cos_pi_f32(a: f32) -> (f32, f32) {
     let c = c.mul_add(r2, 1.0);
     let c = mulsign_f32(c, sx);
 
-    let r3 = r2 * r;
-    let s = S[3];
-    let s = s.mul_add(r2, S[2]);
+    let s = S[2];
     let s = s.mul_add(r2, S[1]);
-    let s = r.mul_add(S[0], r3 * s);
+    let s = s.mul_add(r2, S[0]);
+    let s = r.mul_add(std::f32::consts::PI, r * r2.mul_add(s, -8.742278e-8));
 
     let (s, c) = if i & 1 != 0 { (c, s) } else { (s, c) };
 
diff --git a/narcissus-maths/tests/exhaustive_f32.rs b/narcissus-maths/tests/exhaustive_f32.rs
index ce723a9..3d392f4 100644
--- a/narcissus-maths/tests/exhaustive_f32.rs
+++ b/narcissus-maths/tests/exhaustive_f32.rs
@@ -424,7 +424,7 @@ pub fn exhaustive_sin_pi() {
     let errors = check_exhaustive_f32(ref_sin_pi_f32, |a| sin_cos_pi_f32(a).0, false);
     println!("SIN: {:?}", errors);
     assert_eq!(errors.max_error_ulp, 1);
-    assert_eq!(errors.num_errors, 744_496_052);
+    assert_eq!(errors.num_errors, 55_943_962);
 }
 
 #[test]