1//===-- Implementation header for 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___SUPPORT_MATH_COSF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_COSF16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "sincosf16_utils.h"
17#include "src/__support/FPUtil/FEnvImpl.h"
18#include "src/__support/FPUtil/FPBits.h"
19#include "src/__support/FPUtil/cast.h"
20#include "src/__support/FPUtil/except_value_utils.h"
21#include "src/__support/FPUtil/multiply_add.h"
22#include "src/__support/macros/optimization.h"
23
24namespace LIBC_NAMESPACE_DECL {
25
26namespace math {
27
28LIBC_INLINE constexpr float16 cosf16(float16 x) {
29#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
30 constexpr size_t N_EXCEPTS = 4;
31
32 constexpr fputil::ExceptValues<float16, N_EXCEPTS> COSF16_EXCEPTS{.values: {
33 // (input, RZ output, RU offset, RD offset, RN offset)
34 {.input: 0x2b7c, .rnd_towardzero_result: 0x3bfc, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
35 {.input: 0x4ac1, .rnd_towardzero_result: 0x38b5, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
36 {.input: 0x5c49, .rnd_towardzero_result: 0xb8c6, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 0},
37 {.input: 0x7acc, .rnd_towardzero_result: 0xa474, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 0},
38 }};
39#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
40
41 using namespace sincosf16_internal;
42 using FPBits = fputil::FPBits<float16>;
43 FPBits xbits(x);
44
45 uint16_t x_u = xbits.uintval();
46 uint16_t x_abs = x_u & 0x7fff;
47 float xf = x;
48
49 // Range reduction:
50 // For |x| > pi/32, we perform range reduction as follows:
51 // Find k and y such that:
52 // x = (k + y) * pi/32
53 // k is an integer, |y| < 0.5
54 //
55 // This is done by performing:
56 // k = round(x * 32/pi)
57 // y = x * 32/pi - k
58 //
59 // Once k and y are computed, we then deduce the answer by the cosine of sum
60 // formula:
61 // cos(x) = cos((k + y) * pi/32)
62 // = cos(k * pi/32) * cos(y * pi/32) -
63 // sin(k * pi/32) * sin(y * pi/32)
64
65#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
66 // Handle exceptional values
67 if (auto r = COSF16_EXCEPTS.lookup(x_bits: x_abs); LIBC_UNLIKELY(r.has_value()))
68 return r.value();
69#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
70
71 // cos(+/-0) = 1
72 if (LIBC_UNLIKELY(x_abs == 0U))
73 return fputil::cast<float16>(x: 1.0f);
74
75 // cos(+/-inf) = NaN, and cos(NaN) = NaN
76 if (xbits.is_inf_or_nan()) {
77 if (xbits.is_signaling_nan()) {
78 fputil::raise_except_if_required(FE_INVALID);
79 return FPBits::quiet_nan().get_val();
80 }
81
82 if (xbits.is_inf()) {
83 fputil::set_errno_if_required(EDOM);
84 fputil::raise_except_if_required(FE_INVALID);
85 }
86
87 return x + FPBits::quiet_nan().get_val();
88 }
89
90 float sin_k = 0.0f, cos_k = 0.0f, sin_y = 0.0f, cosm1_y = 0.0f;
91 sincosf16_eval(xf, sin_k, cos_k, sin_y, cosm1_y);
92 // Since, cosm1_y = cos_y - 1, therefore:
93 // cos(x) = cos_k * cos_y - sin_k * sin_y
94 // = cos_k * (cos_y - 1 + 1) - sin_k * sin_y
95 // = cos_k * cosm1_y - sin_k * sin_y + cos_k
96 return fputil::cast<float16>(x: fputil::multiply_add(
97 x: cos_k, y: cosm1_y, z: fputil::multiply_add(x: -sin_k, y: sin_y, z: cos_k)));
98}
99
100} // namespace math
101
102} // namespace LIBC_NAMESPACE_DECL
103
104#endif // LIBC_TYPES_HAS_FLOAT16
105
106#endif // LLVM_LIBC_SRC___SUPPORT_MATH_COSF16_H
107