Skip to content

Commit 0d59443

Browse files
wldfngrswldfngrs
committed
Add cospif16 function
1 parent f5bd36a commit 0d59443

File tree

11 files changed

+296
-0
lines changed

11 files changed

+296
-0
lines changed

‎libc/config/linux/aarch64/entrypoints.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -607,6 +607,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
607607
libc.src.math.canonicalizef16
608608
libc.src.math.ceilf16
609609
libc.src.math.copysignf16
610+
libc.src.math.cospif16
610611
# TODO: aarch64 bug
611612
# Please see https://github.com/llvm/llvm-project/pull/100632#issuecomment-2258772681
612613
# libc.src.math.expf16

‎libc/config/linux/x86_64/entrypoints.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -611,6 +611,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
611611
libc.src.math.ceilf16
612612
libc.src.math.copysignf16
613613
libc.src.math.coshf16
614+
libc.src.math.cospif16
614615
libc.src.math.exp10f16
615616
libc.src.math.exp10m1f16
616617
libc.src.math.exp2f16

‎libc/newhdrgen/yaml/math.yaml

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -206,6 +206,13 @@ functions:
206206
return_type: float
207207
arguments:
208208
- type: float
209+
- name: cospif16
210+
standards:
211+
- stdc
212+
return_type: _Float16
213+
arguments:
214+
- type: _Float16
215+
guard: LIBC_TYPES_HAS_FLOAT16
209216
- name: ddivl
210217
standards:
211218
- stdc

‎libc/src/math/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -95,6 +95,7 @@ add_math_entrypoint_object(coshf)
9595
add_math_entrypoint_object(coshf16)
9696

9797
add_math_entrypoint_object(cospif)
98+
add_math_entrypoint_object(cospif16)
9899

99100
add_math_entrypoint_object(daddl)
100101
add_math_entrypoint_object(daddf128)

‎libc/src/math/cospif16.h

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
//===-- Implementation header for cospif16 ---------------------*- C++ -*-===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
// ===--------------------------------------------------------------------===//
8+
9+
#ifndef LLVM_LIBC_SRC_MATH_COSPIF16_H
10+
#define LLVM_LIBC_SRC_MATH_COSPIF16_H
11+
12+
#include "src/__support/macros/config.h"
13+
#include "src/__support/macros/properties/types.h"
14+
15+
namespace LIBC_NAMESPACE_DECL {
16+
17+
float16 cospif16(float16 x);
18+
19+
} // namespace LIBC_NAMESPACE_DECL
20+
21+
#endif // LLVM_LIBC_SRC_MATH_SINPIF16_H

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

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -422,6 +422,26 @@ add_entrypoint_object(
422422
-O3
423423
)
424424

425+
426+
add_entrypoint_object(
427+
cospif16
428+
SRCS
429+
cospif16.cpp
430+
HDRS
431+
../cospif16.h
432+
DEPENDS
433+
libc.src.__support.common
434+
libc.src.__support.FPUtil.cast
435+
libc.src.__support.FPUtil.fenv_impl
436+
libc.src.__support.FPUtil.fp_bits
437+
libc.src.__support.FPUtil.multiply_add
438+
libc.src.__support.FPUtil.nearest_integer
439+
libc.src.__support.FPUtil.polyeval
440+
libc.src.__support.macros.properties.types
441+
COMPILE_OPTIONS
442+
-O3
443+
)
444+
425445
add_entrypoint_object(
426446
sin
427447
SRCS

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

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
//===-- Half-precision cospif function ------------------------------------===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
//===----------------------------------------------------------------------===//
8+
9+
#include "src/math/cospif16.h"
10+
#include "src/__support/FPUtil/FEnvImpl.h"
11+
#include "src/__support/FPUtil/FPBits.h"
12+
#include "src/__support/FPUtil/PolyEval.h"
13+
#include "src/__support/FPUtil/cast.h"
14+
#include "src/__support/FPUtil/multiply_add.h"
15+
#include "src/__support/FPUtil/nearest_integer.h"
16+
#include "src/__support/common.h"
17+
#include "src/__support/macros/config.h"
18+
19+
namespace LIBC_NAMESPACE_DECL {
20+
// Lookup table for sin(k * pi / 32) with k = 0, ..., 63.
21+
// Table is generated with Sollya as follows:
22+
// > display = hexadecimal;
23+
// > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); };
24+
static constexpr float SIN_K_PI_OVER_32[64] = {
25+
0x0.0p0, 0x1.917a6cp-4, 0x1.8f8b84p-3, 0x1.294062p-2,
26+
0x1.87de2ap-2, 0x1.e2b5d4p-2, 0x1.1c73b4p-1, 0x1.44cf32p-1,
27+
0x1.6a09e6p-1, 0x1.8bc806p-1, 0x1.a9b662p-1, 0x1.c38b3p-1,
28+
0x1.d906bcp-1, 0x1.e9f416p-1, 0x1.f6297cp-1, 0x1.fd88dap-1,
29+
0x1p0, 0x1.fd88dap-1, 0x1.f6297cp-1, 0x1.e9f416p-1,
30+
0x1.d906bcp-1, 0x1.c38b3p-1, 0x1.a9b662p-1, 0x1.8bc806p-1,
31+
0x1.6a09e6p-1, 0x1.44cf32p-1, 0x1.1c73b4p-1, 0x1.e2b5d4p-2,
32+
0x1.87de2ap-2, 0x1.294062p-2, 0x1.8f8b84p-3, 0x1.917a6cp-4,
33+
0x0.0p0, -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294062p-2,
34+
-0x1.87de2ap-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, -0x1.44cf32p-1,
35+
-0x1.6a09e6p-1, -0x1.8bc806p-1, -0x1.a9b662p-1, -0x1.c38b3p-1,
36+
-0x1.d906bcp-1, -0x1.e9f416p-1, -0x1.f6297ep-1, -0x1.fd88dap-1,
37+
-0x1p0, -0x1.fd88dap-1, -0x1.f6297cp-1, -0x1.e9f416p-1,
38+
-0x1.d906bcp-1, -0x1.c38b3p-1, -0x1.a9b662p-1, -0x1.8bc806p-1,
39+
-0x1.6a09e6p-1, -0x1.44cf32p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2,
40+
-0x1.87de2ap-2, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4};
41+
42+
static LIBC_INLINE int32_t range_reduction(float x, float &y) {
43+
float kf = fputil::nearest_integer(x * 32);
44+
y = fputil::multiply_add<float>(x, 32.0, -kf);
45+
46+
return static_cast<int32_t>(kf);
47+
}
48+
49+
LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) {
50+
using FPBits = typename fputil::FPBits<float16>;
51+
FPBits xbits(x);
52+
53+
uint16_t x_u = xbits.uintval();
54+
uint16_t x_abs = x_u & 0x7fff;
55+
56+
// Range reduction:
57+
// For |x| > 1/32, we perform range reduction as follows:
58+
// Find k and y such that:
59+
// x = (k + y) * 1/32
60+
// k is an integer
61+
// |y| < 0.5
62+
//
63+
// This is done by performing:
64+
// k = round(x * 32)
65+
// y = x * 32 - k
66+
//
67+
// Once k and y are computed, we then deduce the answer by the sine of sum
68+
// formula:
69+
// sin(x * pi) = sin((k + y) * pi/32)
70+
// = sin(k * pi/32) * cos(y * pi/32) + sin (y * pi/32) * cos (k *
71+
// pi/32)
72+
// The values of sin(k * pi/32) and cos (k * pi/32) for k = 0...63 are
73+
// precomputed and stored using a vector of 64 single precision floats. sin(y
74+
// * pi/32) and cos(y * pi/32) are computed using degree-9 chebyshev
75+
// polynomials generated by Sollya.
76+
77+
// For signed zeros
78+
if (LIBC_UNLIKELY(x_abs == 0U)) return fputil::cast<float16>(1.0f);
79+
80+
// Numbers greater or equal to 2^10 are integers, or infinity, or NaN
81+
if (LIBC_UNLIKELY(x_abs >= 0x6400)) {
82+
if (LIBC_UNLIKELY(x_abs <= 0x67FF)) {
83+
return fputil::cast<float16>((x_abs & 0x1) ? -1.0f : 1.0f);
84+
}
85+
86+
// Check for NaN or infintiy values
87+
if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
88+
// If value is equal to infinity
89+
if (x_abs == 0x7c00) {
90+
fputil::set_errno_if_required(EDOM);
91+
fputil::raise_except_if_required(FE_INVALID);
92+
}
93+
94+
return x + FPBits::quiet_nan().get_val();
95+
}
96+
97+
return fputil::cast<float16>(1.0f);
98+
}
99+
100+
101+
float f32 = x;
102+
float y;
103+
int32_t k = range_reduction(f32, y);
104+
105+
float sin_k = SIN_K_PI_OVER_32[k & 63];
106+
float cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
107+
108+
// Recall;
109+
// cos(x * pi/32) = cos((k + y) * pi/32)
110+
// = cos(y * pi/32) * cos(k * pi/32)
111+
// - sin(y * pi/32) * sin(k * pi/32)
112+
// Recall, after range reduction, -0.5 <= y <= 0.5. For very small
113+
// values of y, calculating sin(y * p/32) can be inaccurate. Generating a
114+
// polynomial for sin(y * p/32)/y instead significantly reduces the relative
115+
// errors.
116+
float ysq = y * y;
117+
118+
// Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
119+
// with:
120+
// > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]);
121+
float sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f,
122+
0x1.a03354p-21f, -0x1.ad02d2p-20f);
123+
124+
// Note that cosm1_y = cos(y*pi/32) - 1 = cos_y - 1
125+
// Derivation:
126+
// cos(x * pi) = cos((k + y) * pi/32)
127+
// = cos_k * cos_y + sin_k * sin_y
128+
// = cos_k * (1 + cos_y - 1) + sin_k * sin_y
129+
// Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
130+
// > P = fpminimax(cos(y*pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]);
131+
float cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f,
132+
0x1.a6f7a2p-29f);
133+
134+
if (LIBC_UNLIKELY(sin_y == 0 && cos_k == 0))
135+
return fputil::cast<float16>(0.0f);
136+
137+
// Since, cosm1_y = cos_y - 1, therefore:
138+
// cos(x * pi) = cos_k(cosm1_y) + cos_k - sin_k * sin_y
139+
return fputil::cast<float16>(fputil::multiply_add(cos_k, cosm1_y, fputil::multiply_add(-sin_k, sin_y, cos_k)));
140+
}
141+
} // namespace LIBC_NAMESPACE_DECL

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

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,17 @@ add_fp_unittest(
4545
)
4646

4747

48+
add_fp_unittest(
49+
cospif16_test
50+
NEED_MPFR
51+
SUITE
52+
libc-math-unittests
53+
SRCS
54+
cospif16_test.cpp
55+
DEPENDS
56+
libc.src.math.cospif16
57+
)
58+
4859
add_fp_unittest(
4960
daddl_test
5061
NEED_MPFR

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

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
//===-- Exhaustive test for cospif16 --------------------------------------===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
//===---------------------------------------------------------------------===//
8+
9+
#include "src/math/cospif16.h"
10+
#include "test/UnitTest/FPMatcher.h"
11+
#include "test/UnitTest/Test.h"
12+
#include "utils/MPFRWrapper/MPFRUtils.h"
13+
14+
using LlvmLibcCospif16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
15+
16+
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
17+
18+
// Range: [0, Inf]
19+
static constexpr uint16_t POS_START = 0x0000U;
20+
static constexpr uint16_t POS_STOP = 0x7c00U;
21+
22+
// Range: [-Inf, 0]
23+
static constexpr uint16_t NEG_START = 0x8000U;
24+
static constexpr uint16_t NEG_STOP = 0xfc00U;
25+
26+
TEST_F(LlvmLibcCospif16Test, PositiveRange) {
27+
for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
28+
float16 x = FPBits(v).get_val();
29+
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cospi, x, LIBC_NAMESPACE::cospif16(x), 0.5);
30+
}
31+
}
32+
33+
TEST_F(LlvmLibcCospif16Test, NegativeRange) {
34+
for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
35+
float16 x = FPBits(v).get_val();
36+
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cospi, x, LIBC_NAMESPACE::cospif16(x), 0.5);
37+
}
38+
}

‎libc/test/src/math/smoke/CMakeLists.txt

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,17 @@ add_fp_unittest(
2525
libc.src.__support.FPUtil.fp_bits
2626
)
2727

28+
add_fp_unittest(
29+
cospif16_test
30+
SUITE
31+
libc-math-smoke-tests
32+
SRCS
33+
cospif16_test.cpp
34+
DEPENDS
35+
libc.src.errno.errno
36+
libc.src.math.cospif16
37+
)
38+
2839
add_fp_unittest(
2940
sinf_test
3041
SUITE

0 commit comments

Comments
 (0)