git: 6fe5d4d8c3c7 - stable/13 - Fixes for bugs in sinpi/cospi/tanpi

From: Konstantin Belousov <kib_at_FreeBSD.org>
Date: Thu, 10 Aug 2023 02:58:29 UTC
The branch stable/13 has been updated by kib:

URL: https://cgit.FreeBSD.org/src/commit/?id=6fe5d4d8c3c7cf3a2a81f38300766a733029c763

commit 6fe5d4d8c3c7cf3a2a81f38300766a733029c763
Author:     Steve Kargl <kargl@FreeBSD.org>
AuthorDate: 2023-07-31 22:34:48 +0000
Commit:     Konstantin Belousov <kib@FreeBSD.org>
CommitDate: 2023-08-10 02:57:29 +0000

    Fixes for bugs in sinpi/cospi/tanpi
    
    PR:     272742
    
    (cherry picked from commit 2d3b0a687b910c84606e4bc36176945ad5c60406)
---
 lib/msun/ld128/s_cospil.c   | 50 ++++++++++++++++++++----------------------
 lib/msun/ld128/s_sinpil.c   | 38 ++++++++++++--------------------
 lib/msun/ld128/s_tanpil.c   | 39 ++++++++++++++++-----------------
 lib/msun/ld80/s_cospil.c    | 22 ++++++-------------
 lib/msun/ld80/s_sinpil.c    | 16 +++-----------
 lib/msun/ld80/s_tanpil.c    | 34 ++++++++++++-----------------
 lib/msun/src/math_private.h | 53 +++++++++++++++++++++++++++++++++++++++++++++
 lib/msun/src/s_cospi.c      | 20 +++++------------
 lib/msun/src/s_cospif.c     | 16 ++++++--------
 lib/msun/src/s_sinpi.c      | 13 ++---------
 lib/msun/src/s_sinpif.c     | 10 +++------
 lib/msun/src/s_tanpi.c      | 41 +++++++++++++++++------------------
 lib/msun/src/s_tanpif.c     | 24 ++++++++++----------
 13 files changed, 183 insertions(+), 193 deletions(-)

diff --git a/lib/msun/ld128/s_cospil.c b/lib/msun/ld128/s_cospil.c
index 71acc4485f7b..b21f879c3e84 100644
--- a/lib/msun/ld128/s_cospil.c
+++ b/lib/msun/ld128/s_cospil.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017-2021 Steven G. Kargl
+ * Copyright (c) 2017-2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -28,6 +28,7 @@
  * See ../src/s_cospi.c for implementation details.
  */
 
+#include "fpmath.h"
 #include "math.h"
 #include "math_private.h"
 
@@ -46,8 +47,7 @@ volatile static const double vzero = 0;
 long double
 cospil(long double x)
 {
-	long double ax, c, xf;
-	uint32_t ix;
+	long double ai, ar, ax, c;
 
 	ax = fabsl(x);
 
@@ -72,41 +72,37 @@ cospil(long double x)
 	}
 
 	if (ax < 0x1p112) {
-		/* Split x = n + r with 0 <= r < 1. */
-		xf = (ax + 0x1p112L) - 0x1p112L;        /* Integer part */
-		ax -= xf;                               /* Remainder */
-		if (ax < 0) {
-			ax += 1;
-			xf -= 1;
-		}
+		/* Split ax = ai + ar with 0 <= ar < 1. */
+		FFLOORL128(ax, ai, ar);
 
-		if (ax < 0.5) {
-			if (ax < 0.25)
-				c = ax == 0 ? 1 : __kernel_cospil(ax);
+		if (ar < 0.5) {
+			if (ar < 0.25)
+				c = ar == 0 ? 1 : __kernel_cospil(ar);
 			else
-				c = __kernel_sinpil(0.5 - ax);
+				c = __kernel_sinpil(0.5 - ar);
 		} else {
-			if (ax < 0.75) {
-				if (ax == 0.5)
+			if (ar < 0.75) {
+				if (ar == 0.5)
 					return (0);
-				c = -__kernel_sinpil(ax - 0.5);
+				c = -__kernel_sinpil(ar - 0.5);
 			} else
-				c = -__kernel_cospil(1 - ax);
+				c = -__kernel_cospil(1 - ar);
 		}
-
-		if (xf > 0x1p64)
-			xf -= 0x1p64;
-		if (xf > 0x1p32)
-			xf -= 0x1p32;
-		ix = (uint32_t)xf;
-		return (ix & 1 ? -c : c);
+		return (fmodl(ai, 2.L) == 0 ? c : -c);
 	}
 
 	if (isinf(x) || isnan(x))
 		return (vzero / vzero);
 
 	/*
-	 * |x| >= 0x1p112 is always an even integer, so return 1.
+	 * For |x| >= 0x1p113, it is always an even integer, so return 1.
 	 */
-	return (1);
+	if (ax >= 0x1p113)
+		return (1);
+	/*
+	 * For 0x1p112 <= |x| < 0x1p113 need to determine if x is an even
+	 * or odd integer to return 1 or -1.
+	 */
+
+	return (fmodl(ax, 2.L) == 0 ? 1 : -1);
 }
diff --git a/lib/msun/ld128/s_sinpil.c b/lib/msun/ld128/s_sinpil.c
index cdfa2bcac3ef..c8c205449557 100644
--- a/lib/msun/ld128/s_sinpil.c
+++ b/lib/msun/ld128/s_sinpil.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017-2021 Steven G. Kargl
+ * Copyright (c) 2017-2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -28,6 +28,7 @@
  * See ../src/s_sinpi.c for implementation details.
  */
 
+#include "fpmath.h"
 #include "math.h"
 #include "math_private.h"
 
@@ -46,8 +47,7 @@ volatile static const double vzero = 0;
 long double
 sinpil(long double x)
 {
-	long double ax, hi, lo, s, xf, xhi, xlo;
-	uint32_t ix;
+	long double ai, ar, ax, hi, lo, s, xhi, xlo;
 
 	ax = fabsl(x);
 
@@ -78,35 +78,25 @@ sinpil(long double x)
 	}
 
 	if (ax < 0x1p112) {
-		/* Split x = n + r with 0 <= r < 1. */
-		xf = (ax + 0x1p112L) - 0x1p112L;        /* Integer part */
-		ax -= xf;                               /* Remainder */
-		if (ax < 0) {
-			ax += 1;
-			xf -= 1;
-		}
+		/* Split ax = ai + ar with 0 <= ar < 1. */
+		FFLOORL128(ax, ai, ar);
 
-		if (ax == 0) {
+		if (ar == 0) {
 			s = 0;
 		} else {
-			if (ax < 0.5) {
-				if (ax <= 0.25)
-					s = __kernel_sinpil(ax);
+			if (ar < 0.5) {
+				if (ar <= 0.25)
+					s = __kernel_sinpil(ar);
 				else
-					s = __kernel_cospil(0.5 - ax);
+					s = __kernel_cospil(0.5 - ar);
 			} else {
-				if (ax < 0.75)
-					s = __kernel_cospil(ax - 0.5);
+				if (ar < 0.75)
+					s = __kernel_cospil(ar - 0.5);
 				else
-					s = __kernel_sinpil(1 - ax);
+					s = __kernel_sinpil(1 - ar);
 			}
 
-			if (xf > 0x1p64)
-				xf -= 0x1p64;
-			if (xf > 0x1p32)
-				xf -= 0x1p32;
-			ix = (uint32_t)xf;
-			if (ix & 1) s = -s;
+			s = fmodl(ai, 2.L) == 0 ? s : -s;
 		}
 		return (x < 0 ? -s : s);
 	}
diff --git a/lib/msun/ld128/s_tanpil.c b/lib/msun/ld128/s_tanpil.c
index 90f4aea5c629..2d253bb9f478 100644
--- a/lib/msun/ld128/s_tanpil.c
+++ b/lib/msun/ld128/s_tanpil.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017-2021 Steven G. Kargl
+ * Copyright (c) 2017-2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -28,6 +28,7 @@
  * See ../src/s_tanpi.c for implementation details.
  */
 
+#include "fpmath.h"
 #include "math.h"
 #include "math_private.h"
 
@@ -69,8 +70,8 @@ volatile static const double vzero = 0;
 long double
 tanpil(long double x)
 {
-	long double ax, hi, lo, xf, t;
-	uint32_t ix;
+	long double ai, ar, ax, hi, lo, t;
+	double odd;
 
 	ax = fabsl(x);
 
@@ -88,27 +89,22 @@ tanpil(long double x)
 			}
 			t = __kernel_tanpil(ax);
 		} else if (ax == 0.5)
-			return ((ax - ax) / (ax - ax));
+			t = 1 / vzero;
 		else
 			t = -__kernel_tanpil(1 - ax);
 		return (x < 0 ? -t : t);
 	}
 
-	if (ix < 0x1p112) {
-		/* Split x = n + r with 0 <= r < 1. */
-		xf = (ax + 0x1p112L) - 0x1p112L;        /* Integer part */
-		ax -= xf;                               /* Remainder */
-		if (ax < 0) {
-			ax += 1;
-			xf -= 1;
-		}
-
-		if (ax < 0.5)
-			t = ax == 0 ? 0 : __kernel_tanpil(ax);
-		else if (ax == 0.5)
-			return ((ax - ax) / (ax - ax));
+	if (ax < 0x1p112) {
+		/* Split ax = ai + ar with 0 <= ar < 1. */
+		FFLOORL128(ax, ai, ar);
+		odd = fmodl(ai, 2.L) == 0 ? 1 : -1;
+		if (ar < 0.5)
+			t = ar == 0 ? copysign(0., odd) : __kernel_tanpil(ar);
+		else if (ar == 0.5)
+			t = odd / vzero;
 		else
-			t = -__kernel_tanpil(1 - ax);
+			t = -__kernel_tanpil(1 - ar);
 		return (x < 0 ? -t : t);
 	}
 
@@ -117,7 +113,10 @@ tanpil(long double x)
 		return (vzero / vzero);
 
 	/*
-	 * |x| >= 0x1p112 is always an integer, so return +-0.
+	 * For 0x1p112 <= |x| < 0x1p113 need to determine if x is an even
+	 * or odd integer to set t = +0 or -0.
+	 * For |x| >= 0x1p113, it is always an even integer, so t = 0.
 	 */
-	return (copysignl(0, x));
+	t = fmodl(ax,2.L) == 0  ? 0 : copysign(0., -1.);
+	return (copysignl(t, x));
 }
diff --git a/lib/msun/ld80/s_cospil.c b/lib/msun/ld80/s_cospil.c
index 199479e9eaf9..69620d2f2f33 100644
--- a/lib/msun/ld80/s_cospil.c
+++ b/lib/msun/ld80/s_cospil.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017, 2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -80,18 +80,8 @@ cospil(long double x)
 		RETURNI(c);
 	}
 
-	if (ix < 0x403e) {		/* 1 <= |x| < 0x1p63 */
-		/* Determine integer part of ax. */
-		j0 = ix - 0x3fff + 1;
-		if (j0 < 32) {
-			lx = (lx >> 32) << 32;
-			lx &= ~(((lx << 32)-1) >> j0);
-		} else {
-			m = (uint64_t)-1 >> (j0 + 1);
-			if (lx & m) lx &= ~m;
-		}
-		INSERT_LDBL80_WORDS(x, ix, lx);
-
+	if (ix < 0x403e) {			/* 1 <= |x| < 0x1p63 */
+		FFLOORL80(x, j0, ix, lx);	/* Integer part of ax. */
 		ax -= x;
 		EXTRACT_LDBL80_WORDS(ix, lx, ax);
 
@@ -123,7 +113,9 @@ cospil(long double x)
 		RETURNI(vzero / vzero);
 
 	/*
-	 * |x| >= 0x1p63 is always an even integer, so return 1.
+	 * For 0x1p63 <= |x| < 0x1p64 need to determine if x is an even
+	 * or odd integer to return t = +1 or -1.
+	 * For |x| >= 0x1p64, it is always an even integer, so t = 1.
 	 */
-	RETURNI(1);
+	RETURNI(ix >= 0x403f ? 1 : ((lx & 1) ? -1 : 1));
 }
diff --git a/lib/msun/ld80/s_sinpil.c b/lib/msun/ld80/s_sinpil.c
index 4cefa92352e1..7d9008f9e18f 100644
--- a/lib/msun/ld80/s_sinpil.c
+++ b/lib/msun/ld80/s_sinpil.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017, 2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -88,18 +88,8 @@ sinpil(long double x)
 		RETURNI((hx & 0x8000) ? -s : s);
 	}
 
-	if (ix < 0x403e) {		/* 1 <= |x| < 0x1p63 */
-		/* Determine integer part of ax. */
-		j0 = ix - 0x3fff + 1;
-		if (j0 < 32) {
-			lx = (lx >> 32) << 32;
-			lx &= ~(((lx << 32)-1) >> j0);
-		} else {
-			m = (uint64_t)-1 >> (j0 + 1);
-			if (lx & m) lx &= ~m;
-		}
-		INSERT_LDBL80_WORDS(x, ix, lx);
-
+	if (ix < 0x403e) {			/* 1 <= |x| < 0x1p63 */
+		FFLOORL80(x, j0, ix, lx);	/* Integer part of ax. */
 		ax -= x;
 		EXTRACT_LDBL80_WORDS(ix, lx, ax);
 
diff --git a/lib/msun/ld80/s_tanpil.c b/lib/msun/ld80/s_tanpil.c
index 02451e562025..2d640413af6c 100644
--- a/lib/msun/ld80/s_tanpil.c
+++ b/lib/msun/ld80/s_tanpil.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017, 2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -72,7 +72,7 @@ volatile static const double vzero = 0;
 long double
 tanpil(long double x)
 {
-	long double ax, hi, lo, t;
+	long double ax, hi, lo, odd, t;
 	uint64_t lx, m;
 	uint32_t j0;
 	uint16_t hx, ix;
@@ -98,31 +98,22 @@ tanpil(long double x)
 			}
 			t = __kernel_tanpil(ax);
 		} else if (ax == 0.5)
-			RETURNI((ax - ax) / (ax - ax));
+			t = 1 / vzero;
 		else
 			t = -__kernel_tanpil(1 - ax);
 		RETURNI((hx & 0x8000) ? -t : t);
 	}
 
-	if (ix < 0x403e) {		/* 1 <= |x| < 0x1p63 */
-		/* Determine integer part of ax. */
-		j0 = ix - 0x3fff + 1;
-		if (j0 < 32) {
-			lx = (lx >> 32) << 32;
-			lx &= ~(((lx << 32)-1) >> j0);
-		} else {
-			m = (uint64_t)-1 >> (j0 + 1);
-			if (lx & m) lx &= ~m;
-		}
-		INSERT_LDBL80_WORDS(x, ix, lx);
-
+	if (ix < 0x403e) {			/* 1 <= |x| < 0x1p63 */
+		FFLOORL80(x, j0, ix, lx);	/* Integer part of ax. */
+		odd = (uint64_t)x & 1 ? -1 : 1;
 		ax -= x;
 		EXTRACT_LDBL80_WORDS(ix, lx, ax);
 
 		if (ix < 0x3ffe)		/* |x| < 0.5 */
-			t = ax == 0 ? 0 : __kernel_tanpil(ax);
-		else if (ax == 0.5)
-			RETURNI((ax - ax) / (ax - ax));
+			t = ix == 0 ? copysignl(0, odd) : __kernel_tanpil(ax);
+		else if (ax == 0.5L)
+			t = odd / vzero;
 		else
 			t = -__kernel_tanpil(1 - ax);
 		RETURNI((hx & 0x8000) ? -t : t);
@@ -133,7 +124,10 @@ tanpil(long double x)
 		RETURNI(vzero / vzero);
 
 	/*
-	 * |x| >= 0x1p63 is always an integer, so return +-0.
+	 * For 0x1p63 <= |x| < 0x1p64 need to determine if x is an even
+	 * or odd integer to set t = +0 or -0.
+	 * For |x| >= 0x1p64, it is always an even integer, so t = 0.
 	 */
-	RETURNI(copysignl(0, x));
+	t = ix >= 0x403f ? 0 : (copysignl(0, (lx & 1) ? -1 : 1));
+	RETURNI((hx & 0x8000) ? -t : t);
 }
diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h
index 429b9c70b215..814a05ebce67 100644
--- a/lib/msun/src/math_private.h
+++ b/lib/msun/src/math_private.h
@@ -688,6 +688,59 @@ irintl(long double x)
 }
 #endif
 
+/*
+ * The following are fast floor macros for 0 <= |x| < 0x1p(N-1), where
+ * N is the precision of the type of x. These macros are used in the
+ * half-cycle trignometric functions (e.g., sinpi(x)).
+ */
+#define	FFLOORF(x, j0, ix) do {			\
+	(j0) = (((ix) >> 23) & 0xff) - 0x7f;	\
+	(ix) &= ~(0x007fffff >> (j0));		\
+	SET_FLOAT_WORD((x), (ix));		\
+} while (0)
+
+#define	FFLOOR(x, j0, ix, lx) do {				\
+	(j0) = (((ix) >> 20) & 0x7ff) - 0x3ff;			\
+	if ((j0) < 20) {					\
+		(ix) &= ~(0x000fffff >> (j0));			\
+		(lx) = 0;					\
+	} else {						\
+		(lx) &= ~((uint32_t)0xffffffff >> ((j0) - 20));	\
+	}							\
+	INSERT_WORDS((x), (ix), (lx));				\
+} while (0)
+
+#define	FFLOORL80(x, j0, ix, lx) do {			\
+	j0 = ix - 0x3fff + 1;				\
+	if ((j0) < 32) {				\
+		(lx) = ((lx) >> 32) << 32;		\
+		(lx) &= ~((((lx) << 32)-1) >> (j0));	\
+	} else {					\
+		uint64_t _m;				\
+		_m = (uint64_t)-1 >> (j0);		\
+		if ((lx) & _m) (lx) &= ~_m;		\
+	}						\
+	INSERT_LDBL80_WORDS((x), (ix), (lx));		\
+} while (0)
+
+#define FFLOORL128(x, ai, ar) do {			\
+	union IEEEl2bits u;				\
+	uint64_t m;					\
+	int e;						\
+	u.e = (x);					\
+	e = u.bits.exp - 16383;				\
+	if (e < 48) {					\
+		m = ((1llu << 49) - 1) >> (e + 1);	\
+		u.bits.manh &= ~m;			\
+		u.bits.manl = 0;			\
+	} else {					\
+		m = (uint64_t)-1 >> (e - 48);		\
+		u.bits.manl &= ~m;			\
+	}						\
+	(ai) = u.e;					\
+	(ar) = (x) - (ai);				\
+} while (0)
+
 #ifdef DEBUG
 #if defined(__amd64__) || defined(__i386__)
 #define	breakpoint()	asm("int $3")
diff --git a/lib/msun/src/s_cospi.c b/lib/msun/src/s_cospi.c
index 2e2f92733a86..f97570dc8792 100644
--- a/lib/msun/src/s_cospi.c
+++ b/lib/msun/src/s_cospi.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017, 2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -104,20 +104,10 @@ cospi(double x)
 	}
 
 	if (ix < 0x43300000) {		/* 1 <= |x| < 0x1p52 */
-		/* Determine integer part of ax. */
-		j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
-		if (j0 < 20) {
-			ix &= ~(0x000fffff >> j0);
-			lx = 0;
-		} else {
-			lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
-		}
-		INSERT_WORDS(x, ix, lx);
-
+		FFLOOR(x, j0, ix, lx);	/* Integer part of ax. */
 		ax -= x;
 		EXTRACT_WORDS(ix, lx, ax);
 
-
 		if (ix < 0x3fe00000) {		/* |x| < 0.5 */
 			if (ix < 0x3fd00000)	/* |x| < 0.25 */
 				c = ix == 0 ? 1 : __kernel_cospi(ax);
@@ -143,9 +133,11 @@ cospi(double x)
 		return (vzero / vzero);
 
 	/*
-	 * |x| >= 0x1p52 is always an even integer, so return 1.
+	 * For 0x1p52 <= |x| < 0x1p53 need to determine if x is an even
+	 * or odd integer to return +1 or -1.
+	 * For |x| >= 0x1p53, it is always an even integer, so return 1.
 	 */
-	return (1);
+	return (ix < 0x43400000 ? ((lx & 1) ? -1 : 1) : 1);
 }
 
 #if LDBL_MANT_DIG == 53
diff --git a/lib/msun/src/s_cospif.c b/lib/msun/src/s_cospif.c
index 4dd881395baf..44d19f165025 100644
--- a/lib/msun/src/s_cospif.c
+++ b/lib/msun/src/s_cospif.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017,2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -71,12 +71,8 @@ cospif(float x)
 		return (c);
 	}
 
-	if (ix < 0x4b000000) {			/* 1 <= |x| < 0x1p23 */
-		/* Determine integer part of ax. */
-		j0 = ((ix >> 23) & 0xff) - 0x7f;
-		ix &= ~(0x007fffff >> j0);
-		SET_FLOAT_WORD(x, ix);
-
+	if (ix < 0x4b000000) {		/* 1 <= |x| < 0x1p23 */
+		FFLOORF(x, j0, ix);	/* Integer part of ax. */
 		ax -= x;
 		GET_FLOAT_WORD(ix, ax);
 
@@ -103,7 +99,9 @@ cospif(float x)
 		return (vzero / vzero);
 
 	/*
-	 * |x| >= 0x1p23 is always an even integer, so return 1.
+	 * For 0x1p23 <= |x| < 0x1p24 need to determine if x is an even
+	 * or odd integer to return +1 or -1.
+	 * For |x| >= 0x1p24, it is always an even integer, so return 1.
 	 */
-	return (1);
+	return (ix < 0x4b800000 ? ((ix & 1) ? -1 : 1) : 1);
 }
diff --git a/lib/msun/src/s_sinpi.c b/lib/msun/src/s_sinpi.c
index bc3759e567a3..8b388de863c3 100644
--- a/lib/msun/src/s_sinpi.c
+++ b/lib/msun/src/s_sinpi.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017, 2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -118,16 +118,7 @@ sinpi(double x)
 	}
 
 	if (ix < 0x43300000) {			/* 1 <= |x| < 0x1p52 */
-		/* Determine integer part of ax. */
-		j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
-		if (j0 < 20) {
-			ix &= ~(0x000fffff >> j0);
-			lx = 0;
-		} else {
-			lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
-		}
-		INSERT_WORDS(x, ix, lx);
-
+		FFLOOR(x, j0, ix, lx);	/* Integer part of ax. */
 		ax -= x;
 		EXTRACT_WORDS(ix, lx, ax);
 
diff --git a/lib/msun/src/s_sinpif.c b/lib/msun/src/s_sinpif.c
index c9f76f8a2358..21082dee7d9c 100644
--- a/lib/msun/src/s_sinpif.c
+++ b/lib/msun/src/s_sinpif.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017,2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -81,12 +81,8 @@ sinpif(float x)
 		return ((hx & 0x80000000) ? -s : s);
 	}
 
-	if (ix < 0x4b000000) {			/* 1 <= |x| < 0x1p23 */
-		/* Determine integer part of ax. */
-		j0 = ((ix >> 23) & 0xff) - 0x7f;
-		ix &= ~(0x007fffff >> j0);
-		SET_FLOAT_WORD(x, ix);
-
+	if (ix < 0x4b000000) {		/* 1 <= |x| < 0x1p23 */
+		FFLOORF(x, j0, ix);	/* Integer part of ax. */
 		ax -= x;
 		GET_FLOAT_WORD(ix, ax);
 
diff --git a/lib/msun/src/s_tanpi.c b/lib/msun/src/s_tanpi.c
index f911d56156b3..cd00adbcb86e 100644
--- a/lib/msun/src/s_tanpi.c
+++ b/lib/msun/src/s_tanpi.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017, 2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -56,11 +56,15 @@
  * 5. Special cases:
  *
  *    tanpi(+-0) = +-0
- *    tanpi(+-n) = +-0, for positive integers n.
+ *    tanpi(n) = +0 for positive even and negative odd integer n.
+ *    tanpi(n) = -0 for positive odd and negative even integer n.
  *    tanpi(+-n+1/4) = +-1, for positive integers n.
- *    tanpi(+-n+1/2) = NaN, for positive integers n.
- *    tanpi(+-inf) = NaN.  Raises the "invalid" floating-point exception.
- *    tanpi(nan) = NaN.  Raises the "invalid" floating-point exception.
+ *    tanpi(n+1/2) = +inf and raises the FE_DIVBYZERO exception for 
+ *    even integers n.   
+ *    tanpi(n+1/2) = -inf and raises the FE_DIVBYZERO exception for 
+ *    odd integers n.   
+ *    tanpi(+-inf) = NaN and raises the FE_INVALID exception.
+ *    tanpi(nan) = NaN and raises the FE_INVALID exception.
  */
 
 #include <float.h>
@@ -106,7 +110,7 @@ volatile static const double vzero = 0;
 double
 tanpi(double x)
 {
-	double ax, hi, lo, t;
+	double ax, hi, lo, odd, t;
 	uint32_t hx, ix, j0, lx;
 
 	EXTRACT_WORDS(hx, lx, x);
@@ -132,30 +136,22 @@ tanpi(double x)
 			}
 			t = __kernel_tanpi(ax);
 		} else if (ax == 0.5)
-			return ((ax - ax) / (ax - ax));
+			t = 1 / vzero;
 		else
 			t = - __kernel_tanpi(1 - ax);
 		return ((hx & 0x80000000) ? -t : t);
 	}
 
 	if (ix < 0x43300000) {		/* 1 <= |x| < 0x1p52 */
-		/* Determine integer part of ax. */
-		j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
-		if (j0 < 20) {
-			ix &= ~(0x000fffff >> j0);
-			lx = 0;
-		} else {
-			lx &= ~(((uint32_t)(0xffffffff)) >> (j0 - 20));
-		}
-		INSERT_WORDS(x,ix,lx);
-
+		FFLOOR(x, j0, ix, lx);	/* Integer part of ax. */
+		odd = (uint64_t)x & 1 ? -1 : 1;
 		ax -= x;
 		EXTRACT_WORDS(ix, lx, ax);
 
 		if (ix < 0x3fe00000)		/* |x| < 0.5 */
-			t = ax == 0 ? 0 : __kernel_tanpi(ax);
+			t = ix == 0 ? copysign(0, odd) : __kernel_tanpi(ax);
 		else if (ax == 0.5)
-			return ((ax - ax) / (ax - ax));
+			t = odd / vzero;
 		else
 			t = - __kernel_tanpi(1 - ax);
 
@@ -167,9 +163,12 @@ tanpi(double x)
 		return (vzero / vzero);
 
 	/*
-	 * |x| >= 0x1p52 is always an integer, so return +-0.
+	 * For 0x1p52 <= |x| < 0x1p53 need to determine if x is an even
+	 * or odd integer to set t = +0 or -0.
+	 * For |x| >= 0x1p54, it is always an even integer, so t = 0.
 	 */
-	return (copysign(0, x));
+	t = ix >= 0x43400000 ? 0 : (copysign(0, (lx & 1) ? -1 : 1));
+	return ((hx & 0x80000000) ? -t : t);
 }
 
 #if LDBL_MANT_DIG == 53
diff --git a/lib/msun/src/s_tanpif.c b/lib/msun/src/s_tanpif.c
index 6d4b627d1cf9..12dd8f838976 100644
--- a/lib/msun/src/s_tanpif.c
+++ b/lib/msun/src/s_tanpif.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017,2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -58,7 +58,7 @@ volatile static const float vzero = 0;
 float
 tanpif(float x)
 {
-	float ax, hi, lo, t;
+	float ax, hi, lo, odd, t;
 	uint32_t hx, ix, j0;
 
 	GET_FLOAT_WORD(hx, x);
@@ -79,25 +79,22 @@ tanpif(float x)
 			}
 			t = __kernel_tanpif(ax);
 		} else if (ix == 0x3f000000)
-			return ((ax - ax) / (ax - ax));
+			t = 1 / vzero;
 		else
 			t = - __kernel_tanpif(1 - ax);
 		return ((hx & 0x80000000) ? -t : t);
 	}
 
 	if (ix < 0x4b000000) {		/* 1 <= |x| < 0x1p23 */
-		/* Determine integer part of ax. */
-		j0 = ((ix >> 23) & 0xff) - 0x7f;
-		ix &= ~(0x007fffff >> j0);
-		SET_FLOAT_WORD(x, ix);
-
+		FFLOORF(x, j0, ix);	/* Integer part of ax. */
+		odd = (uint32_t)x & 1 ? -1 : 1;
 		ax -= x;
 		GET_FLOAT_WORD(ix, ax);
 
 		if (ix < 0x3f000000)		/* |x| < 0.5 */
-			t = ix == 0 ? 0 : __kernel_tanpif(ax);
+			t = ix == 0 ? copysignf(0, odd) : __kernel_tanpif(ax);
 		else if (ix == 0x3f000000)
-			return ((ax - ax) / (ax - ax));
+			t = odd / vzero;
 		else
 			t = - __kernel_tanpif(1 - ax);
 		return ((hx & 0x80000000) ? -t : t);
@@ -108,7 +105,10 @@ tanpif(float x)
 		return (vzero / vzero);
 
 	/*
-	 * |x| >= 0x1p23 is always an integer, so return +-0.
+	 * For 0x1p23 <= |x| < 0x1p24 need to determine if x is an even
+	 * or odd integer to set t = +0 or -0.
+	 * For |x| >= 0x1p24, it is always an even integer, so t = 0.
 	 */
-	return (copysignf(0, x));
+	t = ix >= 0x4b800000 ? 0 : (copysignf(0, (ix & 1) ? -1 : 1));
+	return ((hx & 0x80000000) ? -t : t);
 }