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

[via libc-commits]

================
@@ -0,0 +1,84 @@
+//===-- Collection of utils for sinf16/cosf16 -------------------*- C++ -*-===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+// Lookup table for sin(k * pi / 32) with k = 0, ..., 63.
+// Table is generated with Sollya as follows:
+// > display = hexadecimmal;
+// > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); };
+constexpr float SIN_K_PI_OVER_32[64] = {
+    0x0.0p0,        0x1.917a6cp-4,  0x1.8f8b84p-3,  0x1.294062p-2,
+    0x1.87de2ap-2,  0x1.e2b5d4p-2,  0x1.1c73b4p-1,  0x1.44cf32p-1,
+    0x1.6a09e6p-1,  0x1.8bc806p-1,  0x1.a9b662p-1,  0x1.c38b3p-1,
+    0x1.d906bcp-1,  0x1.e9f416p-1,  0x1.f6297cp-1,  0x1.fd88dap-1,
+    0x1p0,          0x1.fd88dap-1,  0x1.f6297cp-1,  0x1.e9f416p-1,
+    0x1.d906bcp-1,  0x1.c38b3p-1,   0x1.a9b662p-1,  0x1.8bc806p-1,
+    0x1.6a09e6p-1,  0x1.44cf32p-1,  0x1.1c73b4p-1,  0x1.e2b5d4p-2,
+    0x1.87de2ap-2,  0x1.294062p-2,  0x1.8f8b84p-3,  0x1.917a6cp-4,
+    0x0.0p0,        -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294062p-2,
+    -0x1.87de2ap-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, -0x1.44cf32p-1,
+    -0x1.6a09e6p-1, -0x1.8bc806p-1, -0x1.a9b662p-1, -0x1.c38b3p-1,
+    -0x1.d906bcp-1, -0x1.e9f416p-1, -0x1.f6297ep-1, -0x1.fd88dap-1,
+    -0x1p0,         -0x1.fd88dap-1, -0x1.f6297cp-1, -0x1.e9f416p-1,
+    -0x1.d906bcp-1, -0x1.c38b3p-1,  -0x1.a9b662p-1, -0x1.8bc806p-1,
+    -0x1.6a09e6p-1, -0x1.44cf32p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2,
+    -0x1.87de2ap-2, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4};
+
+LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) {
+  float kf = fputil::nearest_integer(x * 32);
+  y = fputil::multiply_add<float>(x, 32.0, -kf);
+
+  return static_cast<int32_t>(kf);
+}
+
+LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k,
+                                  float &sin_y, float &cosm1_y) {
+  float y;
+  int32_t k = range_reduction_sincospif16(xf, y);
+
+  sin_k = SIN_K_PI_OVER_32[k & 63];
+  cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
+
+  // Recall;
+  // sin(x * pi/32) = sin((k + y) * pi/32)
+  // 		    = sin(y * pi/32) * cos(k * pi/32) + cos(y * pi/32) * sin(k *
+  // pi/32) Recall, after range reduction, -0.5 <= y <= 0.5. For very small
+  // values of y, calculating sin(y * p/32) can be inaccurate. Generating a
+  // polynomial for sin(y * p/32)/y instead significantly reduces the relative
+  // errors.
----------------
overmighty wrote:

Nit: formatting.

```suggestion
  // sin(x * pi/32) = sin((k + y) * pi/32)
  //                = sin(y * pi/32) * cos(k * pi/32) +
  //                  cos(y * pi/32) * sin(k * pi/32)
  // Recall, after range reduction, -0.5 <= y <= 0.5. For very small values of
  // y, calculating sin(y * p/32) can be inaccurate. Generating a polynomial for
  // sin(y * p/32)/y instead significantly reduces the relative errors.
```

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