Skip to content

Commit 1d5f809

Browse files
wldfngrswldfngrs
committed
Fixed rounding errors
1 parent 49b868b commit 1d5f809

File tree

6 files changed

+57
-29
lines changed

6 files changed

+57
-29
lines changed

‎libc/src/math/generic/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -509,6 +509,7 @@ add_entrypoint_object(
509509
libc.src.__support.FPUtil.cast
510510
libc.src.__support.FPUtil.fenv_impl
511511
libc.src.__support.FPUtil.fp_bits
512+
libc.src.__support.FPUtil.except_value_utils
512513
libc.src.__support.FPUtil.multiply_add
513514
libc.src.__support.macros.optimization
514515
COMPILE_OPTIONS

‎libc/src/math/generic/sincosf16_utils.h‎

Lines changed: 10 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@
1111

1212
#include "src/__support/FPUtil/FPBits.h"
1313
#include "src/__support/FPUtil/PolyEval.h"
14+
#include "src/__support/FPUtil/cast.h"
1415
#include "src/__support/FPUtil/nearest_integer.h"
1516
#include "src/__support/common.h"
1617
#include "src/__support/macros/config.h"
@@ -39,11 +40,6 @@ constexpr float SIN_K_PI_OVER_32[64] = {
3940
-0x1.6a09e6p-1, -0x1.44cf32p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2,
4041
-0x1.87de2ap-2, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4};
4142

42-
static constexpr float THIRTYTWO_OVER_PI[4] = {
43-
//0x1.45fp3, 0x1.837p-12, -0x1.b1b8p-28, -0x1.f568p-43
44-
0x1.46p3, -0x1.9f4p-10, 0x1.b94p-22, -0x1.bcp-36
45-
};
46-
4743
LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) {
4844
float kf = fputil::nearest_integer(x * 32);
4945
y = fputil::multiply_add<float>(x, 32.0, -kf);
@@ -52,17 +48,17 @@ LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) {
5248
}
5349

5450
LIBC_INLINE int32_t range_reduction_sincosf16(float x, float &y) {
55-
double xd = x;
56-
double prod = xd * 0x1.45f306dc9c883p3;
51+
double prod = x * 0x1.45f306dc9c883p3;
5752
double kf = fputil::nearest_integer(prod);
58-
double yd = prod - kf;
53+
y = static_cast<float>(prod - kf);
5954

60-
y = static_cast<float>(yd);
6155
return static_cast<int32_t>(kf);
6256
}
6357

64-
static LIBC_INLINE void sincosf16_poly_eval(int32_t k, float y, float &sin_k, float &cos_k, float &sin_y, float &cosm1_y) {
65-
58+
static LIBC_INLINE void sincosf16_poly_eval(int32_t k, float y, float &sin_k,
59+
float &cos_k, float &sin_y,
60+
float &cosm1_y) {
61+
6662
sin_k = SIN_K_PI_OVER_32[k & 63];
6763
cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
6864

@@ -84,7 +80,8 @@ static LIBC_INLINE void sincosf16_poly_eval(int32_t k, float y, float &sin_k, fl
8480
0x1.a6f7a2p-29f);
8581
}
8682

87-
LIBC_INLINE void sincosf16_eval(float xf, float &sin_k, float &cos_k, float &sin_y, float &cosm1_y) {
83+
LIBC_INLINE void sincosf16_eval(float xf, float &sin_k, float &cos_k,
84+
float &sin_y, float &cosm1_y) {
8885
float y;
8986
int32_t k = range_reduction_sincosf16(xf, y);
9087

@@ -95,7 +92,7 @@ LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k,
9592
float &sin_y, float &cosm1_y) {
9693
float y;
9794
int32_t k = range_reduction_sincospif16(xf, y);
98-
95+
9996
sincosf16_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);
10097
}
10198

‎libc/src/math/generic/sinf16.cpp‎

Lines changed: 39 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -13,19 +13,28 @@
1313
#include "src/__support/FPUtil/FEnvImpl.h"
1414
#include "src/__support/FPUtil/FPBits.h"
1515
#include "src/__support/FPUtil/cast.h"
16+
#include "src/__support/FPUtil/except_value_utils.h"
1617
#include "src/__support/FPUtil/multiply_add.h"
1718
#include "src/__support/macros/optimization.h"
1819

1920
namespace LIBC_NAMESPACE_DECL {
2021

22+
constexpr size_t N_EXCEPTS = 2;
23+
24+
constexpr fputil::ExceptValues<float16, N_EXCEPTS> SINF16_EXCEPTS{{
25+
// (input, RZ output, RU offset, RD offset, RN offset)
26+
{0x585c, 0x3ba3, 1, 0, 1},
27+
{0x5cb0, 0xbbff, 0, 1, 0},
28+
}};
29+
2130
LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) {
2231
using FPBits = typename fputil::FPBits<float16>;
2332
FPBits xbits(x);
2433

2534
uint16_t x_u = xbits.uintval();
2635
uint16_t x_abs = x_u & 0x7fff;
2736
float xf = x;
28-
37+
2938
// Range reduction:
3039
// For !x| > pi/32, we perform range reduction as follows:
3140
// Find k and y such that:
@@ -41,15 +50,35 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) {
4150
// sin(x) = sin((k + y) * pi/32)
4251
// = sin(k * pi/32) * cos(y * pi/32) +
4352
// sin(y * pi/32) * cos(k * pi/32)
44-
53+
54+
// Handle exceptional values
55+
if (LIBC_UNLIKELY(x_abs == 0x585C || x_abs == 0x5cb0)) {
56+
bool x_sign = x_u >> 15;
57+
if (auto r = SINF16_EXCEPTS.lookup_odd(x_abs, x_sign);
58+
LIBC_UNLIKELY(r.has_value()))
59+
return r.value();
60+
}
61+
62+
int rounding = fputil::quick_get_round();
4563
if (LIBC_UNLIKELY(x_abs <= 0x13d0)) {
46-
int rounding = fputil::quick_get_round();
47-
48-
// For signed zeros
49-
if ((LIBC_UNLIKELY(x_abs == 0U)) ||
50-
(rounding == FE_UPWARD && xbits.is_pos()) ||
51-
(rounding == FE_DOWNWARD && xbits.is_neg()))
64+
if (LIBC_UNLIKELY(x_abs == 0U))
5265
return x;
66+
67+
if ((rounding == FE_UPWARD && xbits.is_pos()) ||
68+
(rounding == FE_DOWNWARD && xbits.is_neg()))
69+
return x;
70+
71+
if (rounding == FE_UPWARD && xbits.is_neg()) {
72+
x_u--;
73+
return FPBits(x_u).get_val();
74+
}
75+
}
76+
77+
if (LIBC_UNLIKELY(x_abs == 0x2b45)) {
78+
if (rounding == FE_DOWNWARD && xbits.is_neg()) {
79+
x_u--;
80+
return FPBits(x_u).get_val();
81+
}
5382
}
5483

5584
if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
@@ -61,7 +90,6 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) {
6190

6291
return x + FPBits::quiet_nan().get_val();
6392
}
64-
6593

6694
float sin_k, cos_k, sin_y, cosm1_y;
6795
sincosf16_eval(xf, sin_k, cos_k, sin_y, cosm1_y);
@@ -71,7 +99,8 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) {
7199

72100
// Since, cosm1_y = cos_y - 1, therfore:
73101
// sin(x) = cos_k * sin_y + sin_k + (cosm1_y * sin_k)
74-
return fputil::cast<float16>(fputil::multiply_add(sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k)));
102+
return fputil::cast<float16>(fputil::multiply_add(
103+
sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k)));
75104
}
76105

77106
} // namespace LIBC_NAMESPACE_DECL

‎libc/src/math/generic/tanpif16.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ namespace LIBC_NAMESPACE_DECL {
2121

2222
constexpr size_t N_EXCEPTS = 21;
2323

24-
constexpr fputil::ExceptValues<float16, N_EXCEPTS> TANF16_EXCEPTS{{
24+
constexpr fputil::ExceptValues<float16, N_EXCEPTS> TANPIF16_EXCEPTS{{
2525
// (input, RZ output, RU offset, RD offset, RN offset)
2626
{0x07f2, 0x0e3d, 1, 0, 0}, {0x086a, 0x0eee, 1, 0, 1},
2727
{0x08db, 0x0fa0, 1, 0, 0}, {0x094c, 0x1029, 1, 0, 0},
@@ -49,7 +49,7 @@ LLVM_LIBC_FUNCTION(float16, tanpif16, (float16 x)) {
4949
return x;
5050

5151
bool x_sign = x_u >> 15;
52-
if (auto r = TANF16_EXCEPTS.lookup_odd(x_abs, x_sign);
52+
if (auto r = TANPIF16_EXCEPTS.lookup_odd(x_abs, x_sign);
5353
LIBC_UNLIKELY(r.has_value()))
5454
return r.value();
5555
}

‎libc/test/src/math/sinf16_test.cpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,13 +26,15 @@ static constexpr uint16_t NEG_STOP = 0xfc00U;
2626
TEST_F(LlvmLibcSinf16Test, PositiveRange) {
2727
for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
2828
float16 x = FPBits(v).get_val();
29-
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, LIBC_NAMESPACE::sinf16(x), 0.5);
29+
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x,
30+
LIBC_NAMESPACE::sinf16(x), 0.5);
3031
}
3132
}
3233

3334
TEST_F(LlvmLibcSinf16Test, NegativeRange) {
3435
for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
3536
float16 x = FPBits(v).get_val();
36-
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, LIBC_NAMESPACE::sinf16(x), 0.5);
37+
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x,
38+
LIBC_NAMESPACE::sinf16(x), 0.5);
3739
}
3840
}

‎libc/test/src/math/smoke/sinf16_test.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,6 @@
66
//
77
//===----------------------------------------------------------------------===//
88

9-
109
#include "src/errno/libc_errno.h"
1110
#include "src/math/sinf16.h"
1211
#include "test/UnitTest/FPMatcher.h"
@@ -16,7 +15,7 @@ using LlvmLibcSinf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
1615

1716
TEST_F(LlvmLibcSinf16Test, SpecialNumbers) {
1817
LIBC_NAMESPACE::libc_errno = 0;
19-
18+
2019
EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinf16(aNaN));
2120
EXPECT_MATH_ERRNO(0);
2221

0 commit comments

Comments
 (0)