Attachment 243812 Details for Bug 272742 – New diff with corrected ld128 support

[patch] New diff with corrected ld128 support

sinpi_new.diff (text/plain), 20.30 KB, created by Steve Kargl on 2023-08-02 22:58:28 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Steve Kargl

Created: 2023-08-02 22:58:28 UTC

Size: 20.30 KB

patch

obsolete

>diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h
>index df526e71e545..33a792de1ae8 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_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..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_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..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));
> }

Actions: View | Diff

Attachments on bug 272742: 243634 | 243812