Skip to content

Commit 697397e

Browse files
wldfngrswldfngrs
committed
added sincosf16_utils header, fixed formatting in existing pr
1 parent bf93bb4 commit 697397e

File tree

6 files changed

+111
-146
lines changed

6 files changed

+111
-146
lines changed

‎libc/src/math/cospif16.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
1-
//===-- Implementation header for cospif16 ---------------------*- C++ -*-===//
1+
//===-- Implementation header for cospif16 ----------------------*- C++ -*-===//
22
//
33
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
44
// See https://llvm.org/LICENSE.txt for license information.
55
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
66
//
7-
// ===--------------------------------------------------------------------===//
7+
//===----------------------------------------------------------------------===//
88

99
#ifndef LLVM_LIBC_SRC_MATH_COSPIF16_H
1010
#define LLVM_LIBC_SRC_MATH_COSPIF16_H

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

Lines changed: 15 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -351,6 +351,18 @@ 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+
.range_reduction
360+
libc.src.__support.FPUtil.fp_bits
361+
libc.src.__support.FPUtil.polyeval
362+
libc.src.__support.FPUtil.nearest_integer
363+
libc.src.__support.common
364+
)
365+
354366
add_header_library(
355367
sincos_eval
356368
HDRS
@@ -422,22 +434,17 @@ add_entrypoint_object(
422434
-O3
423435
)
424436

425-
426437
add_entrypoint_object(
427438
cospif16
428439
SRCS
429440
cospif16.cpp
430441
HDRS
431442
../cospif16.h
432443
DEPENDS
433-
libc.src.__support.common
444+
.sincosf16_utils
434445
libc.src.__support.FPUtil.cast
435446
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
447+
libc.src.__support.FPUtil.fp_bits
441448
COMPILE_OPTIONS
442449
-O3
443450
)
@@ -555,14 +562,10 @@ add_entrypoint_object(
555562
HDRS
556563
../sinpif16.h
557564
DEPENDS
558-
libc.src.__support.common
565+
.sincosf16_utils
559566
libc.src.__support.FPUtil.cast
560567
libc.src.__support.FPUtil.fenv_impl
561568
libc.src.__support.FPUtil.fp_bits
562-
libc.src.__support.FPUtil.multiply_add
563-
libc.src.__support.FPUtil.nearest_integer
564-
libc.src.__support.FPUtil.polyeval
565-
libc.src.__support.macros.properties.types
566569
COMPILE_OPTIONS
567570
-O3
568571
)

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

Lines changed: 4 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -7,51 +7,20 @@
77
//===----------------------------------------------------------------------===//
88

99
#include "src/math/cospif16.h"
10+
#include "sincosf16_utils.h"
1011
#include "src/__support/FPUtil/FEnvImpl.h"
1112
#include "src/__support/FPUtil/FPBits.h"
12-
#include "src/__support/FPUtil/PolyEval.h"
1313
#include "src/__support/FPUtil/cast.h"
1414
#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"
1815

1916
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-
4917
LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) {
5018
using FPBits = typename fputil::FPBits<float16>;
5119
FPBits xbits(x);
5220

5321
uint16_t x_u = xbits.uintval();
5422
uint16_t x_abs = x_u & 0x7fff;
23+
float xf = x;
5524

5625
// Range reduction:
5726
// For |x| > 1/32, we perform range reduction as follows:
@@ -98,38 +67,8 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) {
9867
return fputil::cast<float16>(1.0f);
9968
}
10069

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);
70+
float sin_k, cos_k, sin_y, cosm1_y;
71+
sincospif16_eval(xf, sin_k, cos_k, sin_y, cosm1_y);
13372

13473
if (LIBC_UNLIKELY(sin_y == 0 && cos_k == 0))
13574
return fputil::cast<float16>(0.0f);

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

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
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+
#ifndef LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H
9+
#define LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H
10+
11+
#include "src/__support/FPUtil/FPBits.h"
12+
#include "src/__support/FPUtil/PolyEval.h"
13+
#include "src/__support/FPUtil/nearest_integer.h"
14+
#include "src/__support/common.h"
15+
#include "src/__support/macros/config.h"
16+
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(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, float &sin_y, float &cosm1_y) {
50+
float y;
51+
int32_t k = range_reduction(xf, y);
52+
53+
sin_k = SIN_K_PI_OVER_32[k & 63];
54+
cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
55+
56+
// Recall;
57+
// sin(x * pi/32) = sin((k + y) * pi/32)
58+
// = sin(y * pi/32) * cos(k * pi/32) + cos(y * pi/32) * sin(k *
59+
// pi/32) Recall, after range reduction, -0.5 <= y <= 0.5. For very small
60+
// values of y, calculating sin(y * p/32) can be inaccurate. Generating a
61+
// polynomial for sin(y * p/32)/y instead significantly reduces the relative
62+
// errors.
63+
float ysq = y * y;
64+
65+
// Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
66+
// with: > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]);
67+
sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f,
68+
0x1.a03354p-21f, -0x1.ad02d2p-20f);
69+
70+
// Note that cosm1_y = cos(y*pi/32) - 1 = cos_y - 1
71+
// Derivation:
72+
// sin(x * pi) = sin((k + y) * pi/32)
73+
// = sin_y * cos_k + cos_y * sin_k
74+
// = cos_k * sin_y + sin_k * (1 + cos_y - 1)
75+
// Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
76+
// > P = fpminimax(cos(y*pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]);
77+
cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f,
78+
0x1.a6f7a2p-29f);
79+
}
80+
81+
} // namespace LIBC_NAMESPACE_DECL
82+
83+
#endif // LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H

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

Lines changed: 6 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -7,52 +7,20 @@
77
//===----------------------------------------------------------------------===//
88

99
#include "src/math/sinpif16.h"
10+
#include "sincosf16_utils.h"
1011
#include "src/__support/FPUtil/FEnvImpl.h"
1112
#include "src/__support/FPUtil/FPBits.h"
12-
#include "src/__support/FPUtil/PolyEval.h"
1313
#include "src/__support/FPUtil/cast.h"
1414
#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"
1815

1916
namespace LIBC_NAMESPACE_DECL {
20-
21-
// Lookup table for sin(k * pi / 32) with k = 0, ..., 63.
22-
// Table is generated with Sollya as follows:
23-
// > display = hexadecimmal;
24-
// > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); };
25-
static constexpr float SIN_K_PI_OVER_32[64] = {
26-
0x0.0p0, 0x1.917a6cp-4, 0x1.8f8b84p-3, 0x1.294062p-2,
27-
0x1.87de2ap-2, 0x1.e2b5d4p-2, 0x1.1c73b4p-1, 0x1.44cf32p-1,
28-
0x1.6a09e6p-1, 0x1.8bc806p-1, 0x1.a9b662p-1, 0x1.c38b3p-1,
29-
0x1.d906bcp-1, 0x1.e9f416p-1, 0x1.f6297cp-1, 0x1.fd88dap-1,
30-
0x1p0, 0x1.fd88dap-1, 0x1.f6297cp-1, 0x1.e9f416p-1,
31-
0x1.d906bcp-1, 0x1.c38b3p-1, 0x1.a9b662p-1, 0x1.8bc806p-1,
32-
0x1.6a09e6p-1, 0x1.44cf32p-1, 0x1.1c73b4p-1, 0x1.e2b5d4p-2,
33-
0x1.87de2ap-2, 0x1.294062p-2, 0x1.8f8b84p-3, 0x1.917a6cp-4,
34-
0x0.0p0, -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294062p-2,
35-
-0x1.87de2ap-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, -0x1.44cf32p-1,
36-
-0x1.6a09e6p-1, -0x1.8bc806p-1, -0x1.a9b662p-1, -0x1.c38b3p-1,
37-
-0x1.d906bcp-1, -0x1.e9f416p-1, -0x1.f6297ep-1, -0x1.fd88dap-1,
38-
-0x1p0, -0x1.fd88dap-1, -0x1.f6297cp-1, -0x1.e9f416p-1,
39-
-0x1.d906bcp-1, -0x1.c38b3p-1, -0x1.a9b662p-1, -0x1.8bc806p-1,
40-
-0x1.6a09e6p-1, -0x1.44cf32p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2,
41-
-0x1.87de2ap-2, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4};
42-
43-
static LIBC_INLINE int32_t range_reduction(float x, float &y) {
44-
float kf = fputil::nearest_integer(x * 32);
45-
y = fputil::multiply_add<float>(x, 32.0, -kf);
46-
47-
return static_cast<int32_t>(kf);
48-
}
49-
5017
LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) {
5118
using FPBits = typename fputil::FPBits<float16>;
5219
FPBits xbits(x);
5320

5421
uint16_t x_u = xbits.uintval();
5522
uint16_t x_abs = x_u & 0x7fff;
23+
float xf = x;
5624

5725
// Range reduction:
5826
// For |x| > 1/32, we perform range reduction as follows:
@@ -68,8 +36,8 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) {
6836
// Once k and y are computed, we then deduce the answer by the sine of sum
6937
// formula:
7038
// sin(x * pi) = sin((k + y) * pi/32)
71-
// = sin(k * pi/32) * cos(y * pi/32) + sin (y * pi/32) * cos (k *
72-
// pi/32)
39+
// = sin(k * pi/32) * cos(y * pi/32)
40+
// + sin (y * pi/32) * cos (k * pi/32)
7341
// The values of sin(k * pi/32) and cos (k * pi/32) for k = 0...63 are
7442
// precomputed and stored using a vector of 64 single precision floats. sin(y
7543
// * pi/32) and cos(y * pi/32) are computed using degree-9 chebyshev
@@ -94,36 +62,8 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) {
9462
return FPBits::zero(xbits.sign()).get_val();
9563
}
9664

97-
float f32 = x;
98-
float y;
99-
int32_t k = range_reduction(f32, y);
100-
101-
float sin_k = SIN_K_PI_OVER_32[k & 63];
102-
float cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
103-
104-
// Recall;
105-
// sin(x * pi/32) = sin((k + y) * pi/32)
106-
// = sin(y * pi/32) * cos(k * pi/32) + cos(y * pi/32) * sin(k *
107-
// pi/32) Recall, after range reduction, -0.5 <= y <= 0.5. For very small
108-
// values of y, calculating sin(y * p/32) can be inaccurate. Generating a
109-
// polynomial for sin(y * p/32)/y instead significantly reduces the relative
110-
// errors.
111-
float ysq = y * y;
112-
113-
// Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
114-
// with: > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]);
115-
float sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f,
116-
0x1.a03354p-21f, -0x1.ad02d2p-20f);
117-
118-
// Note that cosm1_y = cos(y*pi/32) - 1 = cos_y - 1
119-
// Derivation:
120-
// sin(x * pi) = sin((k + y) * pi/32)
121-
// = sin_y * cos_k + cos_y * sin_k
122-
// = cos_k * sin_y + sin_k * (1 + cos_y - 1)
123-
// Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
124-
// > P = fpminimax(cos(y*pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]);
125-
float cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f,
126-
0x1.a6f7a2p-29f);
65+
float sin_k, cos_k, sin_y, cosm1_y;
66+
sincospif16_eval(xf, sin_k, cos_k, sin_y, cosm1_y);
12767

12868
if (LIBC_UNLIKELY(sin_y == 0 && sin_k == 0))
12969
return FPBits::zero(xbits.sign()).get_val();

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44
// See https://llvm.org/LICENSE.txt for license information.
55
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
66
//
7-
//===---------------------------------------------------------------------===//
7+
//===----------------------------------------------------------------------===//
88

99
#include "src/math/cospif16.h"
1010
#include "test/UnitTest/FPMatcher.h"

0 commit comments

Comments
 (0)