diff --git a/lib/msun/ld128/s_sinpil.c b/lib/msun/ld128/s_sinpil.c index 39eed9b007b..cdfa2bcac3e 100644 --- a/lib/msun/ld128/s_sinpil.c +++ b/lib/msun/ld128/s_sinpil.c @@ -1,5 +1,5 @@ /*- - * Copyright (c) 2017 Steven G. Kargl + * Copyright (c) 2017-2021 Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -26,9 +26,6 @@ /* * See ../src/s_sinpi.c for implementation details. - * - * FIXME: This has not been compiled nor has it been tested for accuracy. - * FIXME: This should use bit twiddling. */ #include "math.h" @@ -68,7 +65,7 @@ sinpil(long double x) } s = __kernel_sinpil(ax); - return (copysignl(s, x)); + return (x < 0 ? -s : s); } if (ax < 0.5) @@ -77,12 +74,18 @@ sinpil(long double x) s = __kernel_cospil(ax - 0.5); else s = __kernel_sinpil(1 - ax); - return (copysignl(s, x)); + return (x < 0 ? -s : s); } if (ax < 0x1p112) { - xf = floorl(ax); - ax -= xf; + /* 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) { s = 0; } else { @@ -98,14 +101,14 @@ sinpil(long double x) s = __kernel_sinpil(1 - ax); } - if (xf > 0x1p50) - xf -= 0x1p50; - if (xf > 0x1p30) - xf -= 0x1p30; + if (xf > 0x1p64) + xf -= 0x1p64; + if (xf > 0x1p32) + xf -= 0x1p32; ix = (uint32_t)xf; if (ix & 1) s = -s; } - return (copysignl(s, x)); + return (x < 0 ? -s : s); } if (isinf(x) || isnan(x)) diff --git a/lib/msun/ld128/s_cospil.c b/lib/msun/ld128/s_cospil.c index b4bc50bb4d8..71acc4485f7 100644 --- a/lib/msun/ld128/s_cospil.c +++ b/lib/msun/ld128/s_cospil.c @@ -1,5 +1,5 @@ /*- - * Copyright (c) 2017 Steven G. Kargl + * Copyright (c) 2017-2021 Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -26,9 +26,6 @@ /* * See ../src/s_cospi.c for implementation details. - * - * FIXME: This has not been compiled nor has it been tested for accuracy. - * FIXME: This should use bit twiddling. */ #include "math.h" @@ -54,7 +51,7 @@ cospil(long double x) ax = fabsl(x); - if (ax < 1) { + if (ax <= 1) { if (ax < 0.25) { if (ax < 0x1p-60) { if ((int)x == 0) @@ -75,15 +72,21 @@ cospil(long double x) } if (ax < 0x1p112) { - xf = floorl(ax); - ax -= xf; - if (x < 0.5) { - if (x < 0.25) + /* 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) { + if (ax < 0.25) c = ax == 0 ? 1 : __kernel_cospil(ax); else c = __kernel_sinpil(0.5 - ax); } else { - if (x < 0.75) { + if (ax < 0.75) { if (ax == 0.5) return (0); c = -__kernel_sinpil(ax - 0.5); @@ -91,10 +94,10 @@ cospil(long double x) c = -__kernel_cospil(1 - ax); } - if (xf > 0x1p50) - xf -= 0x1p50; - if (xf > 0x1p30) - xf -= 0x1p30; + if (xf > 0x1p64) + xf -= 0x1p64; + if (xf > 0x1p32) + xf -= 0x1p32; ix = (uint32_t)xf; return (ix & 1 ? -c : c); } diff --git a/lib/msun/ld128/s_tanpil.c b/lib/msun/ld128/s_tanpil.c index 7cbbdc8a5d0..90f4aea5c62 100644 --- a/lib/msun/ld128/s_tanpil.c +++ b/lib/msun/ld128/s_tanpil.c @@ -1,5 +1,5 @@ /*- - * Copyright (c) 2017 Steven G. Kargl + * Copyright (c) 2017-2021 Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -26,9 +26,6 @@ /* * See ../src/s_tanpi.c for implementation details. - * - * FIXME: This has not been compiled nor has it been tested for accuracy. - * FIXME: This should use bit twiddling. */ #include "math.h" @@ -75,7 +72,7 @@ tanpil(long double x) long double ax, hi, lo, xf, t; uint32_t ix; - ax = fabsl(ax); + ax = fabsl(x); if (ax < 1) { if (ax < 0.5) { @@ -94,19 +91,25 @@ tanpil(long double x) return ((ax - ax) / (ax - ax)); else t = -__kernel_tanpil(1 - ax); - return (copysignl(t, x)); + return (x < 0 ? -t : t); } if (ix < 0x1p112) { - xf = floorl(ax); - ax -= xf; + /* 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)); else t = -__kernel_tanpil(1 - ax); - return (copysignl(t, x)); + return (x < 0 ? -t : t); } /* x = +-inf or nan. */ @@ -114,7 +117,7 @@ tanpil(long double x) return (vzero / vzero); /* - * |x| >= 0x1p53 is always an integer, so return +-0. + * |x| >= 0x1p112 is always an integer, so return +-0. */ return (copysignl(0, x)); }