Date: Mon, 1 Nov 2021 02:38:54 GMT From: Konstantin Belousov <kib@FreeBSD.org> To: src-committers@FreeBSD.org, dev-commits-src-all@FreeBSD.org, dev-commits-src-main@FreeBSD.org Subject: git: 4f889260c33c - main - sinpi[fl] etc: Fix the ld128 implementations Message-ID: <202111010238.1A12csuF052789@gitrepo.freebsd.org>
next in thread | raw e-mail | index | archive | help
The branch main has been updated by kib: URL: https://cgit.FreeBSD.org/src/commit/?id=4f889260c33c163ab28e0e082b4d7e7562d9c647 commit 4f889260c33c163ab28e0e082b4d7e7562d9c647 Author: Steve Kargl <kargl@FreeBSD.org> AuthorDate: 2021-10-31 22:26:20 +0000 Commit: Konstantin Belousov <kib@FreeBSD.org> CommitDate: 2021-11-01 02:38:19 +0000 sinpi[fl] etc: Fix the ld128 implementations Mark Murray graciously provided access to an aarch64 system to test the ld128 implementations. This patch address * Misuses of copysignl() in sinpil() and tanpil(). * Redo the splitting of argument 'x' into an integer part and remainder. The remainder must satify 0 <= r < 1. * Update the reduction of the integer part to something that can easily be seen as even or odd, e.g., sin(pi*x) = (-1)^n*sin(pi*r) with n <= 2^112 and we an reduce n by subtracting integer powers of 2. * In s_cospil.c, fix typos where 'x' is used where 'ax', the remainder, is required. * In tanpil(), fix the use of an uninitialized variable, ax = fabsl(ax), ax should be x in fabsl(). One item of note, in the limited tested on aarch64, the max ULP for sinpil() and cospil() were less than 1.1 ULP, which is higher that the desired max ULP less than 1. This was traced to the kernel for cosl() in the fundamental interval [0,pi/4]. The coefficients in the minmax polynomial likely need refinement. PR: 218514 MFC after: 1 week --- lib/msun/ld128/s_cospil.c | 31 +++++++++++++++++-------------- lib/msun/ld128/s_sinpil.c | 29 ++++++++++++++++------------- lib/msun/ld128/s_tanpil.c | 23 +++++++++++++---------- 3 files changed, 46 insertions(+), 37 deletions(-) diff --git a/lib/msun/ld128/s_cospil.c b/lib/msun/ld128/s_cospil.c index b4bc50bb4d89..71acc4485f7b 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_sinpil.c b/lib/msun/ld128/s_sinpil.c index 39eed9b007bc..cdfa2bcac3ef 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_tanpil.c b/lib/msun/ld128/s_tanpil.c index 7cbbdc8a5d05..90f4aea5c629 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)); }
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?202111010238.1A12csuF052789>