The branch stable/13 has been updated by kib: URL: https://cgit.FreeBSD.org/src/commit/?id=6fe5d4d8c3c7cf3a2a81f38300766a733029c763
commit 6fe5d4d8c3c7cf3a2a81f38300766a733029c763 Author: Steve Kargl <[email protected]> AuthorDate: 2023-07-31 22:34:48 +0000 Commit: Konstantin Belousov <[email protected]> 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); }
