[libc-commits] [libc] [libc][math] Implement double precision sin correctly rounded to all rounding modes. (PR #95736)
Nick Desaulniers via libc-commits
libc-commits at lists.llvm.org
Tue Jun 18 10:29:25 PDT 2024
================
@@ -0,0 +1,334 @@
+//===-- Range reduction for double precision sin/cos/tan w/ FMA -*- 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_FMA_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_FMA_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 fma {
+
+using fputil::DoubleDouble;
+using Float128 = fputil::DyadicFloat<128>;
+
+LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 32;
+
+// 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};
+LIBC_INLINE constexpr Float128 PI_OVER_128_F128 = {
+ Sign::POS, -133, 0xc90f'daa2'2168'c234'c4c6'628b'80dc'1cd1_u128};
+
+// Digits of 2^(16*i) / pi, generated by Sollya with:
+// For [2..62]:
+// > for i from 3 to 63 do {
+// pi_inv = 2^(16*(i - 3)) / pi;
+// pn = nearestint(pi_inv);
+// pi_frac = pi_inv - pn;
+// a = round(pi_frac, D, RN);
+// b = round(pi_frac - a, D, RN);
+// c = round(pi_frac - a - b, D, RN);
+// d = round(pi_frac - a - b - c, D, RN);
+// print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},");
+// };
+// For [0..1]:
+// The leading bit of 2^(16*(i - 3)) / pi is very small, so we add 0.25 so that
+// the conditions for the algorithms are still satisfied, and one of those
+// conditions guarantees that ulp(0.25 * x_reduced) >= 2, and will safely be
+// discarded.
+// > for i from 0 to 2 do {
+// pi_frac = 0.25 + 2^(16*(i - 3)) / pi;
+// a = round(pi_frac, D, RN);
+// b = round(pi_frac - a, D, RN);
+// c = round(pi_frac - a - b, D, RN);
+// d = round(pi_frac - a - b - c, D, RN);
+// print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},");
+// };
+// For The fast pass using double-double, we only need 3 parts (a, b, c), but
+// for the accurate pass using Float128, instead of using another table of
+// Float128s, we simply add the fourth path (a, b, c, d), which simplify the
+// implementation a bit and saving some memory.
+LIBC_INLINE constexpr double ONE_TWENTY_EIGHT_OVER_PI[64][4] = {
+ {0x1.0000000000014p5, 0x1.7cc1b727220a9p-49, 0x1.3f84eafa3ea6ap-103,
+ -0x1.11f924eb53362p-157},
+ {0x1.0000000145f3p5, 0x1.b727220a94fe1p-49, 0x1.d5f47d4d37703p-104,
+ 0x1.b6295993c439p-158},
+ {0x1.000145f306dcap5, -0x1.bbead603d8a83p-50, 0x1.f534ddc0db629p-106,
+ 0x1.664f10e4107f9p-160},
+ {0x1.45f306dc9c883p5, -0x1.6b01ec5417056p-49, -0x1.6447e493ad4cep-103,
+ 0x1.e21c820ff28b2p-157},
+ {-0x1.f246c6efab581p4, 0x1.3abe8fa9a6eep-53, 0x1.b6c52b3278872p-107,
+ 0x1.07f9458eaf7afp-164},
+ {0x1.391054a7f09d6p4, -0x1.70565911f924fp-53, 0x1.2b3278872084p-107,
+ -0x1.ae9c5421443aap-162},
+ {0x1.529fc2757d1f5p2, 0x1.a6ee06db14acdp-53, -0x1.8778df7c035d4p-107,
+ 0x1.d5ef5de2b0db9p-161},
+ {-0x1.ec54170565912p-1, 0x1.b6c52b3278872p-59, 0x1.07f9458eaf7afp-116,
+ -0x1.d4f246dc8e2dfp-173},
+ {-0x1.505c1596447e5p5, 0x1.b14acc9e21c82p-49, 0x1.fe5163abdebbcp-106,
+ 0x1.586dc91b8e909p-160},
+ {-0x1.596447e493ad5p1, 0x1.93c439041fe51p-54, 0x1.8eaf7aef1586ep-108,
+ -0x1.b7238b7b645a4p-163},
+ {0x1.bb81b6c52b328p5, -0x1.de37df00d74e3p-49, 0x1.7bd778ac36e49p-103,
+ -0x1.1c5bdb22d1ffap-158},
+ {0x1.b6c52b3278872p5, 0x1.07f9458eaf7afp-52, -0x1.d4f246dc8e2dfp-109,
+ 0x1.374b801924bbbp-164},
+ {0x1.2b3278872084p5, -0x1.ae9c5421443aap-50, 0x1.b7246e3a424ddp-106,
+ 0x1.700324977504fp-161},
+ {-0x1.8778df7c035d4p5, 0x1.d5ef5de2b0db9p-49, 0x1.1b8e909374b8p-104,
+ 0x1.924bba8274648p-160},
+ {-0x1.bef806ba71508p4, -0x1.443a9e48db91cp-50, -0x1.6f6c8b47fe6dbp-104,
+ -0x1.115f62e6de302p-158},
+ {-0x1.ae9c5421443aap-2, 0x1.b7246e3a424ddp-58, 0x1.700324977504fp-113,
+ -0x1.cdbc603c429c7p-167},
+ {-0x1.38a84288753c9p5, -0x1.b7238b7b645a4p-51, 0x1.924bba8274648p-112,
+ 0x1.cfe1deb1cb12ap-166},
+ {-0x1.0a21d4f246dc9p3, 0x1.d2126e9700325p-53, -0x1.a22bec5cdbc6p-107,
+ -0x1.e214e34ed658cp-162},
+ {-0x1.d4f246dc8e2dfp3, 0x1.374b801924bbbp-52, -0x1.f62e6de301e21p-106,
+ -0x1.38d3b5963045ep-160},
+ {-0x1.236e4716f6c8bp4, -0x1.1ff9b6d115f63p-50, 0x1.921cfe1deb1cbp-106,
+ 0x1.29a73ee88235fp-162},
+ {0x1.b8e909374b802p4, -0x1.b6d115f62e6dep-50, -0x1.80f10a71a76b3p-105,
+ 0x1.cfba208d7d4bbp-160},
+ {0x1.09374b801924cp4, -0x1.15f62e6de301ep-50, -0x1.0a71a76b2c609p-105,
+ 0x1.1046bea5d7689p-159},
+ {-0x1.68ffcdb688afbp3, -0x1.736f180f10a72p-53, 0x1.62534e7dd1047p-107,
+ -0x1.0568a25dbd8b3p-161},
+ {0x1.924bba8274648p0, 0x1.cfe1deb1cb12ap-54, -0x1.63045df7282b4p-108,
+ -0x1.44bb7b16638fep-162},
+ {-0x1.a22bec5cdbc6p5, -0x1.e214e34ed658cp-50, -0x1.177dca0ad144cp-106,
+ 0x1.213a671c09ad1p-160},
+ {0x1.3a32439fc3bd6p1, 0x1.cb129a73ee882p-54, 0x1.afa975da24275p-109,
+ -0x1.8e3f652e8207p-164},
+ {-0x1.b78c0788538d4p4, 0x1.29a73ee88235fp-50, 0x1.4baed1213a672p-104,
+ -0x1.fb29741037d8dp-159},
+ {0x1.fc3bd63962535p5, -0x1.822efb9415a29p-51, 0x1.a24274ce38136p-105,
+ -0x1.741037d8cdc54p-159},
+ {-0x1.4e34ed658c117p2, -0x1.f7282b4512edfp-52, 0x1.d338e04d68bfp-107,
+ -0x1.bec66e29c67cbp-162},
+ {0x1.62534e7dd1047p5, -0x1.0568a25dbd8b3p-49, -0x1.c7eca5d040df6p-105,
+ -0x1.9b8a719f2b318p-160},
+ {-0x1.63045df7282b4p4, -0x1.44bb7b16638fep-50, 0x1.ad17df904e647p-104,
+ 0x1.639835339f49dp-158},
+ {0x1.d1046bea5d769p5, -0x1.bd8b31c7eca5dp-49, -0x1.037d8cdc538dp-107,
+ 0x1.a99cfa4e422fcp-161},
+ {0x1.afa975da24275p3, -0x1.8e3f652e8207p-52, 0x1.3991d63983534p-106,
+ -0x1.82d8dee81d108p-160},
+ {-0x1.a28976f62cc72p5, 0x1.35a2fbf209cc9p-53, -0x1.4e33e566305b2p-109,
+ 0x1.08bf177bf2507p-163},
+ {-0x1.76f62cc71fb29p5, -0x1.d040df633714ep-49, -0x1.9f2b3182d8defp-104,
+ 0x1.f8bbdf9283b2p-158},
+ {0x1.d338e04d68bfp5, -0x1.bec66e29c67cbp-50, 0x1.9cfa4e422fc5ep-105,
+ -0x1.036be27003b4p-161},
+ {0x1.c09ad17df904ep4, 0x1.91d639835339fp-50, 0x1.272117e2ef7e5p-104,
+ -0x1.7c4e007680022p-158},
+ {0x1.68befc827323bp5, -0x1.c67cacc60b638p-50, 0x1.17e2ef7e4a0ecp-104,
+ 0x1.ff897ffde0598p-158},
+ {-0x1.037d8cdc538dp5, 0x1.a99cfa4e422fcp-49, 0x1.77bf250763ff1p-103,
+ 0x1.7ffde05980fefp-158},
+ {-0x1.8cdc538cf9599p5, 0x1.f49c845f8bbep-50, -0x1.b5f13801da001p-104,
+ 0x1.e05980fef2f12p-158},
+ {-0x1.4e33e566305b2p3, 0x1.08bf177bf2507p-51, 0x1.8ffc4bffef02dp-105,
+ -0x1.fc04343b9d298p-160},
+ {-0x1.f2b3182d8dee8p4, -0x1.d1081b5f13802p-52, 0x1.2fffbc0b301fep-107,
+ -0x1.a1dce94beb25cp-163},
+ {-0x1.8c16c6f740e88p5, -0x1.036be27003b4p-49, -0x1.0fd33f8086877p-109,
+ -0x1.d297d64b824b2p-164},
+ {0x1.3908bf177bf25p5, 0x1.d8ffc4bffef03p-53, -0x1.9fc04343b9d29p-108,
+ -0x1.f592e092c9813p-162},
+ {0x1.7e2ef7e4a0ec8p4, -0x1.da00087e99fcp-56, -0x1.0d0ee74a5f593p-110,
+ 0x1.f6d367ecf27cbp-166},
+ {-0x1.081b5f13801dap4, -0x1.0fd33f8086877p-61, -0x1.d297d64b824b2p-116,
+ -0x1.8130d834f648bp-170},
+ {-0x1.af89c00ed0004p5, -0x1.fa67f010d0ee7p-50, -0x1.297d64b824b26p-104,
+ -0x1.30d834f648b0cp-162},
+ {-0x1.c00ed00043f4dp5, 0x1.fde5e2316b415p-55, -0x1.2e092c98130d8p-110,
+ -0x1.a7b24585ce04dp-165},
+ {0x1.2fffbc0b301fep5, -0x1.a1dce94beb25cp-51, -0x1.25930261b069fp-107,
+ 0x1.b74f463f669e6p-162},
+ {-0x1.0fd33f8086877p3, -0x1.d297d64b824b2p-52, -0x1.8130d834f648bp-106,
+ -0x1.738132c3402bap-163},
+ {-0x1.9fc04343b9d29p4, -0x1.f592e092c9813p-50, -0x1.b069ec9161738p-107,
+ -0x1.32c3402ba515bp-163},
+ {-0x1.0d0ee74a5f593p2, 0x1.f6d367ecf27cbp-54, 0x1.36e9e8c7ecd3dp-111,
+ -0x1.00ae9456c229cp-165},
+ {-0x1.dce94beb25c12p5, -0x1.64c0986c1a7b2p-49, -0x1.161738132c34p-103,
+ -0x1.5d28ad8453814p-158},
+ {-0x1.4beb25c12593p5, -0x1.30d834f648b0cp-50, 0x1.8fd9a797fa8b6p-104,
+ -0x1.5b08a7028341dp-159},
+ {0x1.b47db4d9fb3cap4, -0x1.a7b24585ce04dp-53, 0x1.3cbfd45aea4f7p-107,
+ 0x1.63f5f2f8bd9e8p-161},
+ {-0x1.25930261b069fp5, 0x1.b74f463f669e6p-50, -0x1.5d28ad8453814p-110,
+ -0x1.a0e84c2f8c608p-166},
+ {0x1.fb3c9f2c26dd4p4, -0x1.738132c3402bap-51, -0x1.456c229c0a0dp-105,
+ -0x1.d0985f18c10ebp-159},
+ {-0x1.b069ec9161738p5, -0x1.32c3402ba515bp-51, -0x1.14e050683a131p-108,
+ 0x1.0739f78a5292fp-162},
+ {-0x1.ec9161738132cp5, -0x1.a015d28ad8454p-50, 0x1.faf97c5ecf41dp-104,
+ -0x1.821d6b5b4565p-160},
+ {-0x1.61738132c3403p5, 0x1.16ba93dd63f5fp-49, 0x1.7c5ecf41ce7dep-104,
+ 0x1.4a525d4d7f6bfp-159},
+ {0x1.fb34f2ff516bbp3, -0x1.b08a7028341d1p-51, 0x1.9e839cfbc5295p-105,
+ -0x1.a2b2809409dc1p-159},
+ {0x1.3cbfd45aea4f7p5, 0x1.63f5f2f8bd9e8p-49, 0x1.ce7de294a4baap-104,
+ -0x1.404a04ee072a3p-158},
+ {-0x1.5d28ad8453814p2, -0x1.a0e84c2f8c608p-54, -0x1.d6b5b45650128p-108,
+ -0x1.3b81ca8bdea7fp-164},
+ {-0x1.15b08a7028342p5, 0x1.7b3d0739f78a5p-50, 0x1.497535fdafd89p-105,
+ -0x1.ca8bdea7f33eep-164},
+};
+
+LIBC_INLINE int range_reduction_small(double x, DoubleDouble &u) {
+ double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI[3][0];
+ double kd = fputil::nearest_integer(prod_hi);
+ int k = static_cast<int>(static_cast<int64_t>(kd));
+
+ // Let y = x - k * (pi/128)
+ // Then |y| < pi / 256
+ // With extra rounding errors, we can bound |y| < 2^-6.
+ double y_hi = fputil::multiply_add(kd, -PI_OVER_128.hi, x); // Exact
+ // u_hi + u_lo ~ (y_hi + kd*(-PI_OVER_128[1]))
+ // and |u_lo| < 2* ulp(u_hi)
+ // The upper bound 2^-6 is over-estimated, we should still have:
+ // |u_hi + u_lo| < 2^-6.
+ u.hi = fputil::multiply_add(kd, -PI_OVER_128.lo, y_hi);
+ u.lo = y_hi - u.hi; // Exact;
+ u.lo = fputil::multiply_add(kd, -PI_OVER_128.lo, u.lo);
+
+ return k;
+}
+
+// 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.
+//
+// Note: this algorithm works correctly without FMA instruction for the default
+// rounding mode, round-to-nearest. The limitation is due to Veltkamp's
+// Splitting algorithm used by exact_mult: double x double -> double-double.
+LIBC_INLINE int range_reduction_large(double x, DoubleDouble &u) {
+ // |x| >= 2^32.
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ 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 - 2))
+ xbits.set_biased_exponent((x_e_m62 & 15) + FPBits::EXP_BIAS + 62);
+ double x_reduced = xbits.get_val();
+ // x * c_hi = ph.hi + ph.lo exactly.
+ DoubleDouble ph =
+ fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]);
+ // x * c_mid = pm.hi + pm.lo exactly.
+ DoubleDouble pm =
+ fputil::exact_mult(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);
+ u = fputil::quick_mult(y, PI_OVER_128);
+ int k = static_cast<int>(kh) + static_cast<int>(km);
+
+ return k;
+}
+
+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);
+ Float128 u = fputil::quick_mul(y, PI_OVER_128_F128);
+
+ return u;
+}
+
+// Maybe not redo-ing most of the computation, instead getting
+// y_hi, idx, pm.lo, x_reduced from range_reduction_large.
+LIBC_INLINE Float128 range_reduction_large_f128(double x) {
+ // |x| >= 2^32.
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ 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 - 2))
+ xbits.set_biased_exponent((x_e_m62 & 15) + FPBits::EXP_BIAS + 62);
+ double x_reduced = xbits.get_val();
+ // x * c_hi = ph.hi + ph.lo exactly.
+ DoubleDouble ph =
+ fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]);
+ // x * c_mid = pm.hi + pm.lo exactly.
+ DoubleDouble pm =
+ fputil::exact_mult(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
+ Float128 y_lo_0(x_reduced * ONE_TWENTY_EIGHT_OVER_PI[idx][3]);
+ Float128 y_lo_1 = fputil::quick_mul(
+ Float128(x_reduced), Float128(ONE_TWENTY_EIGHT_OVER_PI[idx][2]));
+ Float128 y_lo_2(pm.lo);
+ Float128 y_hi_f128(y_hi);
+
+ using fputil::quick_add;
+ Float128 y =
+ quick_add(y_hi_f128, quick_add(y_lo_2, quick_add(y_lo_1, y_lo_0)));
+ Float128 u = fputil::quick_mul(y, PI_OVER_128_F128);
+
+ return u;
----------------
nickdesaulniers wrote:
```suggestion
return fputil::quick_mul(y, PI_OVER_128_F128);
```
https://github.com/llvm/llvm-project/pull/95736
More information about the libc-commits
mailing list