[libc-commits] [libc] [libc][math] Implement double precision sin correctly rounded to all rounding modes. (PR #95736)

[Nick Desaulniers via libc-commits]

================
@@ -0,0 +1,67 @@
+//===-- Range reduction for double precision sin/cos/tan --------*- 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_RANGE_REDUCTION_DOUBLE_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_H
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+
+namespace LIBC_NAMESPACE {
+
+using fputil::DoubleDouble;
+
+LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 23;
+
+// Digits of pi/128, generated by Sollya with:
+// > a = round(pi/128, D, RN);
+// > b = round(pi/128 - a, D, RN);
+LIBC_INLINE constexpr DoubleDouble PI_OVER_128 = {0x1.1a62633145c07p-60,
+                                                  0x1.921fb54442d18p-6};
+
+// Digits of -pi/128, generated by Sollya with:
+// > a = round(pi/128, 25, RN);
+// > b = round(pi/128 - a, 23, RN);
+// > c = round(pi/128 - a - b, 25, RN);
+// > d = round(pi/128 - a - b - c, D, RN);
+// The precisions of the parts are chosen so that:
+// 1)  k * a, k * b, k * c are exact in double precision
+// 2)  k * b + fractional part of (k * a) is exact in double precsion
+LIBC_INLINE constexpr double MPI_OVER_128[4] = {
+    -0x1.921fb5p-6, -0x1.110b48p-32, +0x1.ee59dap-56, -0x1.98a2e03707345p-83};
+
+LIBC_INLINE constexpr double ONE_TWENTY_EIGHT_OVER_PI_D = 0x1.45f306dc9c883p5;
+
+namespace generic {
+
+LIBC_INLINE int range_reduction_small(double x, DoubleDouble &u) {
+  double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI_D;
+  double kd = fputil::nearest_integer(prod_hi);
+  int k = static_cast<int>(kd);
+
+  // x - k * (pi/128)
----------------
nickdesaulniers wrote:

Mind expanding this comment to use the terms `c`, `h_hi`, and `y_lo`?

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