Attachment 243634 Details for Bug 272742 – patch to fix half-cycle trigonometric functions

[patch] patch to fix half-cycle trigonometric functions

z.diff (text/plain), 17.80 KB, created by Steve Kargl on 2023-07-26 22:18:32 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Steve Kargl

Created: 2023-07-26 22:18:32 UTC

Size: 17.80 KB

patch

obsolete

>diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h
>index df526e71e545..04f80886ee71 100644
>--- a/lib/msun/src/math_private.h
>+++ b/lib/msun/src/math_private.h
>@@ -688,6 +688,51 @@ 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)
>+
>+/* Split ax into an integer part xf and a remainder ax. */
>+#define	FFLOORL128(ax, xf) do {			\
>+	(xf) = ((ax) + 0x1p112L) - 0x1p112L;	\
>+	(ax) -= (xf);				\
>+	if ((ax) < 0) {				\
>+		(ax) += 1;			\
>+		(xf) -= 1;			\
>+	}					\
>+} while (0)
>+
> #ifdef DEBUG
> #if defined(__amd64__) || defined(__i386__)
> #define	breakpoint()	asm("int $3")
>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_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_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);
> }
>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_sinpi.c b/lib/msun/src/s_sinpi.c
>index bc3759e567a3..392b18aaa0e0 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,15 +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_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/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/ld128/s_cospil.c b/lib/msun/ld128/s_cospil.c
>index 71acc4485f7b..998cc4ba12d4 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
>@@ -72,13 +72,8 @@ 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 = n + r with 0 <= r < 1. */
>+		FFLOORL128(ax, xf);
> 
> 		if (ax < 0.5) {
> 			if (ax < 0.25)
>@@ -106,7 +101,14 @@ cospil(long double 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.
>+	 */
>+	if (xf >= 0x1p113)
>+		return (1);
>+	/*
>+	 * For 0x1p112 <= |x| < 0x1p113 need to determine if x is an even
>+	 * or odd integer to return +1 or -1.
> 	 */
>-	return (1);
>+
>+	return (fmodl(xf, 2.L) == 1 ? -1 : 1);
> }
>diff --git a/lib/msun/ld128/s_sinpil.c b/lib/msun/ld128/s_sinpil.c
>index cdfa2bcac3ef..751dd46cc183 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
>@@ -79,12 +79,7 @@ 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;
>-		}
>+		FFLOORL128(ax, xf);
> 
> 		if (ax == 0) {
> 			s = 0;
>diff --git a/lib/msun/ld128/s_tanpil.c b/lib/msun/ld128/s_tanpil.c
>index 90f4aea5c629..b9ec3ef50f80 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
>@@ -69,7 +69,7 @@ volatile static const double vzero = 0;
> long double
> tanpil(long double x)
> {
>-	long double ax, hi, lo, xf, t;
>+	long double ax, hi, lo, odd, t,  xf;
> 	uint32_t ix;
> 
> 	ax = fabsl(x);
>@@ -88,7 +88,7 @@ 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);
>@@ -96,17 +96,13 @@ tanpil(long double x)
> 
> 	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;
>-		}
>+		FFLOORL128(ax, xf);
>+		odd = fmodl(xf, 2.L) == 1 ? -1 : 1;
> 
> 		if (ax < 0.5)
> 			t = ax == 0 ? 0 : __kernel_tanpil(ax);
> 		else if (ax == 0.5)
>-			return ((ax - ax) / (ax - ax));
>+			t = odd / vzero;
> 		else
> 			t = -__kernel_tanpil(1 - ax);
> 		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(xf, 2.L) == 1 ? 0 : (copysign(0, (lx & 1) ? -1 : 1));
>+	return ((hx & 0x80000000) ? -t : t);
> }

Actions: View | Diff

Attachments on bug 272742: 243634 | 243812