From 63cbb9862bbaedf428ca2701f60ddf7e65f21787 Mon Sep 17 00:00:00 2001
From: gnezdo <gnezdo@openbsd.org>
Date: Fri, 21 Jan 2022 03:12:43 +0000
Subject: [PATCH] Copy cos(3) software implementation from FreeBSD-13

The result passes more tests from msun suite. In particular,
	testacc(cos, -0.34061437849088045332L, 0.94254960031831729956L,
		ALL_STD_EXCEPT, FE_INEXACT);
matches instead of being 1e-16 off.

OK mbuhl@, tested by mbuhl@ and phessler@
---
 lib/libm/src/k_cos.c | 45 ++++++++++++++++++--------------------------
 lib/libm/src/s_cos.c |  6 +++++-
 2 files changed, 23 insertions(+), 28 deletions(-)

diff --git a/lib/libm/src/k_cos.c b/lib/libm/src/k_cos.c
index 8f3882b6a00..0839243e90c 100644
--- a/lib/libm/src/k_cos.c
+++ b/lib/libm/src/k_cos.c
@@ -36,13 +36,17 @@
  *			  ~ cos(x) - x*y,
  *	   a correction term is necessary in cos(x) and hence
  *		cos(x+y) = 1 - (x*x/2 - (r - x*y))
- *	   For better accuracy when x > 0.3, let qx = |x|/4 with
- *	   the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
- *	   Then
- *		cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
- *	   Note that 1-qx and (x*x/2-qx) is EXACT here, and the
- *	   magnitude of the latter is at least a quarter of x*x/2,
- *	   thus, reducing the rounding error in the subtraction.
+ *	   For better accuracy, rearrange to
+ *		cos(x+y) ~ w + (tmp + (r-x*y))
+ *	   where w = 1 - x*x/2 and tmp is a tiny correction term
+ *	   (1 - x*x/2 == w + tmp exactly in infinite precision).
+ *	   The exactness of w + tmp in infinite precision depends on w
+ *	   and tmp having the same precision as x.  If they have extra
+ *	   precision due to compiler bugs, then the extra precision is
+ *	   only good provided it is retained in all terms of the final
+ *	   expression for cos().  Retention happens in all cases tested
+ *	   under FreeBSD, so don't pessimize things by forcibly clipping
+ *	   any extra precision in w.
  */
 
 #include "math.h"
@@ -60,25 +64,12 @@ C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
 double
 __kernel_cos(double x, double y)
 {
-	double a,hz,z,r,qx;
-	int32_t ix;
-	GET_HIGH_WORD(ix,x);
-	ix &= 0x7fffffff;			/* ix = |x|'s high word*/
-	if(ix<0x3e400000) {			/* if x < 2**27 */
-	    if(((int)x)==0) return one;		/* generate inexact */
-	}
+	double hz,z,r,w;
+
 	z  = x*x;
-	r  = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
-	if(ix < 0x3FD33333) 			/* if |x| < 0.3 */ 
-	    return one - (0.5*z - (z*r - x*y));
-	else {
-	    if(ix > 0x3fe90000) {		/* x > 0.78125 */
-		qx = 0.28125;
-	    } else {
-	        INSERT_WORDS(qx,ix-0x00200000,0);	/* x/4 */
-	    }
-	    hz = 0.5*z-qx;
-	    a  = one-qx;
-	    return a - (hz - (z*r-x*y));
-	}
+	w  = z*z;
+	r  = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
+	hz = 0.5*z;
+	w  = one-hz;
+	return w + (((one-w)-hz) + (z*r-x*y));
 }
diff --git a/lib/libm/src/s_cos.c b/lib/libm/src/s_cos.c
index 8b923d5fe61..1406504e9ab 100644
--- a/lib/libm/src/s_cos.c
+++ b/lib/libm/src/s_cos.c
@@ -57,7 +57,11 @@ cos(double x)
 
     /* |x| ~< pi/4 */
 	ix &= 0x7fffffff;
-	if(ix <= 0x3fe921fb) return __kernel_cos(x,z);
+	if(ix <= 0x3fe921fb) {
+	    if(ix<0x3e46a09e)			/* if x < 2**-27 * sqrt(2) */
+		if(((int)x)==0) return 1.0;	/* generate inexact */
+	    return __kernel_cos(x,z);
+	}
 
     /* cos(Inf or NaN) is NaN */
 	else if (ix>=0x7ff00000) return x-x;
-- 
2.20.1