[libc-commits] [libc] [libc][math][c23] Add cospif16 function (PR #113001)

[via libc-commits]

================
@@ -0,0 +1,81 @@
+//===-- Half-precision cospif function ------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/cospif16.h"
+#include "sincosf16_utils.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/multiply_add.h"
+
+namespace LIBC_NAMESPACE_DECL {
+LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) {
+  using FPBits = typename fputil::FPBits<float16>;
+  FPBits xbits(x);
+
+  uint16_t x_u = xbits.uintval();
+  uint16_t x_abs = x_u & 0x7fff;
+  float xf = x;
+
+  // Range reduction:
+  // For |x| > 1/32, we perform range reduction as follows:
+  // Find k and y such that:
+  //   x = (k + y) * 1/32
+  //   k is an integer
+  //   |y| < 0.5
+  //
+  // This is done by performing:
+  //   k = round(x * 32)
+  //   y = x * 32 - k
+  //
+  // Once k and y are computed, we then deduce the answer by the sine of sum
+  // formula:
+  //   cos(x * pi) = cos((k + y) * pi/32)
+  //           = cos(k * pi/32) * cos(y * pi/32)
+  //           + sin(y * pi/32) * sin(k * pi/32)
+  // The values of sin(k * pi/32) and cos (k * pi/32) for k = 0...63 are
+  // precomputed and stored using a vector of 64 single precision floats. sin(y
+  // * pi/32) and cos(y * pi/32) are computed using degree-9 chebyshev
+  // polynomials generated by Sollya.
+
+  // For signed zeros
+  if (LIBC_UNLIKELY(x_abs == 0U))
+    return fputil::cast<float16>(1.0f);
+
+  // Numbers greater or equal to 2^10 are integers, or infinity, or NaN
+  if (LIBC_UNLIKELY(x_abs >= 0x6400)) {
+    if (LIBC_UNLIKELY(x_abs <= 0x67FF)) {
+      return fputil::cast<float16>((x_abs & 0x1) ? -1.0f : 1.0f);
+    }
+
+    // Check for NaN or infintiy values
+    if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
+      // If value is equal to infinity
+      if (x_abs == 0x7c00) {
+        fputil::set_errno_if_required(EDOM);
+        fputil::raise_except_if_required(FE_INVALID);
----------------
overmighty wrote:

Should include "hdr/errno_macros.h" and "hdr/fenv_macros.h" for `EDOM` and `FE_INVALID`.

https://github.com/llvm/llvm-project/pull/113001