git: 6fe5d4d8c3c7 - stable/13 - Fixes for bugs in sinpi/cospi/tanpi
- Go to: [ bottom of page ] [ top of archives ] [ this month ]
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);
}