[libc-commits] [libc] [libc][math] Implement double precision sin correctly rounded to all rounding modes. (PR #95736)
via libc-commits
libc-commits at lists.llvm.org
Fri Jun 21 21:12:15 PDT 2024
================
@@ -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.
----------------
lntue wrote:
I refactored large range reduction functions to a common class to reuse all those 4 variables.
https://github.com/llvm/llvm-project/pull/95736
More information about the libc-commits
mailing list