git: 5033d0b56f5b - stable/13 - sinpi[fl] etc: Fix the ld128 implementations
- Go to: [ bottom of page ] [ top of archives ] [ this month ]
Date: Tue, 02 Nov 2021 01:17:03 UTC
The branch stable/13 has been updated by kib:
URL: https://cgit.FreeBSD.org/src/commit/?id=5033d0b56f5be90fcca4ede92f28f9335eb20b4a
commit 5033d0b56f5be90fcca4ede92f28f9335eb20b4a
Author: Steve Kargl <kargl@FreeBSD.org>
AuthorDate: 2021-10-31 22:26:20 +0000
Commit: Konstantin Belousov <kib@FreeBSD.org>
CommitDate: 2021-11-02 01:16:29 +0000
sinpi[fl] etc: Fix the ld128 implementations
PR: 218514
(cherry picked from commit 4f889260c33c163ab28e0e082b4d7e7562d9c647)
---
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));
}