[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,173 @@
+//===-- 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_COMMON_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_COMMON_H
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+#include "src/__support/integer_literals.h"
+
+namespace LIBC_NAMESPACE {
+
+namespace generic {
+
+using LIBC_NAMESPACE::fputil::DoubleDouble;
+using Float128 = LIBC_NAMESPACE::fputil::DyadicFloat<128>;
+
+LIBC_INLINE constexpr Float128 PI_OVER_128_F128 = {
+    Sign::POS, -133, 0xc90f'daa2'2168'c234'c4c6'628b'80dc'1cd1_u128};
+
+// Note: The look-up tables ONE_TWENTY_EIGHT_OVER_PI is selected to be either
+// from fma:: or nofma:: namespace.
+
+// For large range |x| >= 2^32, we use the exponent of x to find 3 double-chunks
+// of 128/pi c_hi, c_mid, c_lo such that:
+//   1) ulp(round(x * c_hi, D, RN)) >= 256,
+//   2) If x * c_hi = ph_hi + ph_lo and x * c_mid = pm_hi + pm_lo, then
+//        min(ulp(ph_lo), ulp(pm_hi)) >= 2^-53.
+//   3) ulp(round(x * c_lo, D, RN)) <= 2^-7x.
+// This will allow us to do quick computations as:
+//   (x * 256/pi) ~ x * (c_hi + c_mid + c_lo)    (mod 256)
+//                ~ ph_lo + pm_hi + pm_lo + (x * c_lo)
+// Then,
+//   round(x * 128/pi) = round(ph_lo + pm_hi)    (mod 256)
+// And the high part of fractional part of (x * 128/pi) can simply be:
+//   {x * 128/pi}_hi = {ph_lo + pm_hi}.
+// To prevent overflow when x is very large, we simply scale up
+// (c_hi, c_mid, c_lo) by a fixed power of 2 (based on the index) and scale down
+// x by the same amount.
+//
+
+template <bool NO_FMA>
+LIBC_INLINE unsigned range_reduction_large(double x, DoubleDouble &u) {
+
+  // Digits of pi/128, generated by Sollya with:
+  // > a = round(pi/128, D, RN);
+  // > b = round(pi/128 - a, D, RN);
+  constexpr DoubleDouble PI_OVER_128_DD = {0x1.1a62633145c07p-60,
+                                           0x1.921fb54442d18p-6};
+
+  using FPBits = typename fputil::FPBits<double>;
+  FPBits xbits(x);
+
+  // TODO: The extra exponent gap of 62 below can be reduced a bit for non-FMA
+  // with a more careful analysis, which in turn will reduce the error bound for
+  // non-FMA
+  int x_e_m62 = xbits.get_biased_exponent() - (FPBits::EXP_BIAS + 62);
+  int idx = (x_e_m62 >> 4) + 3;
+  // Scale x down by 2^(-(16 * (idx - 3))
+  xbits.set_biased_exponent((x_e_m62 & 15) + FPBits::EXP_BIAS + 62);
+  // 2^62 <= |x_reduced| < 2^(62 + 16) = 2^78
+  double x_reduced = xbits.get_val();
+  // x * c_hi = ph.hi + ph.lo exactly.
+  DoubleDouble ph =
+      fputil::exact_mult<NO_FMA>(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]);
+  // x * c_mid = pm.hi + pm.lo exactly.
+  DoubleDouble pm =
+      fputil::exact_mult<NO_FMA>(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][1]);
+  // Extract integral parts and fractional parts of (ph.lo + pm.hi).
+  double kh = fputil::nearest_integer(ph.lo);
+  double ph_lo_frac = ph.lo - kh; // Exact
+  double km = fputil::nearest_integer(pm.hi + ph_lo_frac);
+  double pm_hi_frac = pm.hi - km; // Exact
+  // x * 128/pi mod 1 ~ y_hi + y_lo = u.hi + u.lo
+  double y_hi = ph_lo_frac + pm_hi_frac; // Exact
+  // y_lo = x * c_lo + pm.lo
+  double y_lo =
+      fputil::multiply_add(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][2], pm.lo);
+  DoubleDouble y = fputil::exact_add(y_hi, y_lo);
+  // Error bound: with {a} denote the fractional part of a, i.e.:
+  //   {a} = a - round(a)
+  // Then,
+  //   | {x * 128/pi} - (y_hi + y_lo) | <  2 * ulp(x_reduced *
+  //                                         * ONE_TWENTY_EIGHT_OVER_PI[idx][2])
+  // For FMA:
+  //   | {x * 128/pi} - (y_hi + y_lo) | <= 2 * 2^77 * 2^-103 * 2^-52
+  //                                    =  2^-77.
+  //   | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-77.
+  //                                      = 2^-82.
+  // For non-FMA:
+  //   | {x * 128/pi} - (y_hi + y_lo) | <= 2 * 2^77 * 2^-99 * 2^-52
+  //                                    =  2^-73.
+  //   | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-73.
+  //                                      = 2^-78.
+  u = fputil::quick_mult<NO_FMA>(y, PI_OVER_128_DD);
+
+  return static_cast<unsigned>(static_cast<int64_t>(kh) +
+                               static_cast<int64_t>(km));
+}
+
+LIBC_INLINE Float128 range_reduction_small_f128(double x) {
+  double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI[3][0];
+  double kd = fputil::nearest_integer(prod_hi);
+
+  Float128 mk_f128(-kd);
+  Float128 x_f128(x);
+  Float128 p_hi =
+      fputil::quick_mul(x_f128, Float128(ONE_TWENTY_EIGHT_OVER_PI[3][0]));
+  Float128 p_mid =
+      fputil::quick_mul(x_f128, Float128(ONE_TWENTY_EIGHT_OVER_PI[3][1]));
+  Float128 p_lo =
+      fputil::quick_mul(x_f128, Float128(ONE_TWENTY_EIGHT_OVER_PI[3][2]));
+  Float128 s_hi = fputil::quick_add(p_hi, mk_f128);
+  Float128 s_lo = fputil::quick_add(p_mid, p_lo);
+  Float128 y = fputil::quick_add(s_hi, s_lo);
+
+  return fputil::quick_mul(y, PI_OVER_128_F128);
+}
+
+// TODO: Maybe not redo-ing most of the computation, instead getting
+//   y_hi, idx, pm.lo, x_reduced from range_reduction_large.
----------------
nickdesaulniers wrote:

Yeah, this code is pretty dense. I think it would make sense to just DRY it up here in this PR.  So a shared helper that took those 4 as non-const references (in-out params), and `double x` as params?

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