Skip to content

Commit 7395ef5

Browse files
wldfngrsovermighty
wldfngrs
and
overmighty
authored
[libc][math][c23] Add cospif16 function (#113001)
Implementation of `cos` for half precision floating point inputs scaled by pi (i.e., `cospi`), correctly rounded for all rounding modes. --------- Co-authored-by: OverMighty <[email protected]>
1 parent a8398bd commit 7395ef5

File tree

15 files changed

+347
-88
lines changed

15 files changed

+347
-88
lines changed

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

+1
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

+1
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/docs/math/index.rst

+1-1
Original file line numberDiff line numberDiff line change
@@ -280,7 +280,7 @@ Higher Math Functions
280280
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
281281
| cosh | |check| | | | |check| | | 7.12.5.4 | F.10.2.4 |
282282
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
283-
| cospi | |check| | | | | | 7.12.4.12 | F.10.1.12 |
283+
| cospi | |check| | | | |check| | | 7.12.4.12 | F.10.1.12 |
284284
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
285285
| dsqrt | N/A | N/A | |check| | N/A | |check|\* | 7.12.14.6 | F.10.11 |
286286
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+

‎libc/newhdrgen/yaml/math.yaml

+7
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: coshf16
210217
standards:
211218
- stdc

‎libc/src/math/CMakeLists.txt

+1
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

+21
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

+35-5
Original file line numberDiff line numberDiff line change
@@ -351,6 +351,17 @@ add_header_library(
351351
libc.src.__support.common
352352
)
353353

354+
add_header_library(
355+
sincosf16_utils
356+
HDRS
357+
sincosf16_utils.h
358+
DEPENDS
359+
libc.src.__support.FPUtil.fp_bits
360+
libc.src.__support.FPUtil.polyeval
361+
libc.src.__support.FPUtil.nearest_integer
362+
libc.src.__support.common
363+
)
364+
354365
add_header_library(
355366
sincos_eval
356367
HDRS
@@ -422,6 +433,25 @@ add_entrypoint_object(
422433
-O3
423434
)
424435

436+
add_entrypoint_object(
437+
cospif16
438+
SRCS
439+
cospif16.cpp
440+
HDRS
441+
../cospif16.h
442+
DEPENDS
443+
.sincosf16_utils
444+
libc.hdr.errno_macros
445+
libc.hdr.fenv_macros
446+
libc.src.__support.FPUtil.cast
447+
libc.src.__support.FPUtil.fenv_impl
448+
libc.src.__support.FPUtil.fp_bits
449+
libc.src.__support.FPUtil.multiply_add
450+
libc.src.__support.macros.optimization
451+
COMPILE_OPTIONS
452+
-O3
453+
)
454+
425455
add_entrypoint_object(
426456
sin
427457
SRCS
@@ -535,14 +565,14 @@ add_entrypoint_object(
535565
HDRS
536566
../sinpif16.h
537567
DEPENDS
538-
libc.src.__support.common
568+
.sincosf16_utils
569+
libc.hdr.errno_macros
570+
libc.hdr.fenv_macros
539571
libc.src.__support.FPUtil.cast
540572
libc.src.__support.FPUtil.fenv_impl
541-
libc.src.__support.FPUtil.fp_bits
573+
libc.src.__support.FPUtil.fp_bits
542574
libc.src.__support.FPUtil.multiply_add
543-
libc.src.__support.FPUtil.nearest_integer
544-
libc.src.__support.FPUtil.polyeval
545-
libc.src.__support.macros.properties.types
575+
libc.src.__support.macros.optimization
546576
COMPILE_OPTIONS
547577
-O3
548578
)

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

+81
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
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 "hdr/errno_macros.h"
11+
#include "hdr/fenv_macros.h"
12+
#include "sincosf16_utils.h"
13+
#include "src/__support/FPUtil/FEnvImpl.h"
14+
#include "src/__support/FPUtil/FPBits.h"
15+
#include "src/__support/FPUtil/cast.h"
16+
#include "src/__support/FPUtil/multiply_add.h"
17+
#include "src/__support/macros/optimization.h"
18+
19+
namespace LIBC_NAMESPACE_DECL {
20+
21+
LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) {
22+
using FPBits = typename fputil::FPBits<float16>;
23+
FPBits xbits(x);
24+
25+
uint16_t x_u = xbits.uintval();
26+
uint16_t x_abs = x_u & 0x7fff;
27+
float xf = x;
28+
29+
// Range reduction:
30+
// For |x| > 1/32, we perform range reduction as follows:
31+
// Find k and y such that:
32+
// x = (k + y) * 1/32
33+
// k is an integer
34+
// |y| < 0.5
35+
//
36+
// This is done by performing:
37+
// k = round(x * 32)
38+
// y = x * 32 - k
39+
//
40+
// Once k and y are computed, we then deduce the answer by the sine of sum
41+
// formula:
42+
// cos(x * pi) = cos((k + y) * pi/32)
43+
// = cos(k * pi/32) * cos(y * pi/32) +
44+
// sin(y * pi/32) * sin(k * pi/32)
45+
46+
// For signed zeros
47+
if (LIBC_UNLIKELY(x_abs == 0U))
48+
return fputil::cast<float16>(1.0f);
49+
50+
// Numbers greater or equal to 2^10 are integers, or infinity, or NaN
51+
if (LIBC_UNLIKELY(x_abs >= 0x6400)) {
52+
if (LIBC_UNLIKELY(x_abs <= 0x67FF))
53+
return fputil::cast<float16>((x_abs & 0x1) ? -1.0f : 1.0f);
54+
55+
// Check for NaN or infintiy values
56+
if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
57+
// If value is equal to infinity
58+
if (x_abs == 0x7c00) {
59+
fputil::set_errno_if_required(EDOM);
60+
fputil::raise_except_if_required(FE_INVALID);
61+
}
62+
63+
return x + FPBits::quiet_nan().get_val();
64+
}
65+
66+
return fputil::cast<float16>(1.0f);
67+
}
68+
69+
float sin_k, cos_k, sin_y, cosm1_y;
70+
sincospif16_eval(xf, sin_k, cos_k, sin_y, cosm1_y);
71+
72+
if (LIBC_UNLIKELY(sin_y == 0 && cos_k == 0))
73+
return fputil::cast<float16>(0.0f);
74+
75+
// Since, cosm1_y = cos_y - 1, therefore:
76+
// cos(x * pi) = cos_k(cosm1_y) + cos_k - sin_k * sin_y
77+
return fputil::cast<float16>(fputil::multiply_add(
78+
cos_k, cosm1_y, fputil::multiply_add(-sin_k, sin_y, cos_k)));
79+
}
80+
81+
} // namespace LIBC_NAMESPACE_DECL

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

+77
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
//===-- Collection of utils for sinf16/cosf16 -------------------*- 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_GENERIC_SINCOSF16_UTILS_H
10+
#define LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H
11+
12+
#include "src/__support/FPUtil/FPBits.h"
13+
#include "src/__support/FPUtil/PolyEval.h"
14+
#include "src/__support/FPUtil/nearest_integer.h"
15+
#include "src/__support/common.h"
16+
#include "src/__support/macros/config.h"
17+
18+
namespace LIBC_NAMESPACE_DECL {
19+
20+
// Lookup table for sin(k * pi / 32) with k = 0, ..., 63.
21+
// Table is generated with Sollya as follows:
22+
// > display = hexadecimmal;
23+
// > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); };
24+
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+
LIBC_INLINE int32_t range_reduction_sincospif16(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+
LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k,
50+
float &sin_y, float &cosm1_y) {
51+
float y;
52+
int32_t k = range_reduction_sincospif16(xf, y);
53+
54+
sin_k = SIN_K_PI_OVER_32[k & 63];
55+
cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
56+
57+
// Recall, after range reduction, -0.5 <= y <= 0.5. For very small values of
58+
// y, calculating sin(y * p/32) can be inaccurate. Generating a polynomial for
59+
// sin(y * p/32)/y instead significantly reduces the relative errors.
60+
float ysq = y * y;
61+
62+
// Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
63+
// with:
64+
// > Q = fpminimax(sin(y * pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]);
65+
sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f,
66+
0x1.a03354p-21f, -0x1.ad02d2p-20f);
67+
68+
// Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya
69+
// with:
70+
// > P = fpminimax(cos(y * pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]);
71+
cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f,
72+
0x1.a6f7a2p-29f);
73+
}
74+
75+
} // namespace LIBC_NAMESPACE_DECL
76+
77+
#endif // LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H

0 commit comments

Comments
 (0)