Skip to content

[libc][math] Improve performance of double precision trig functions. #111793

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 7 commits into from
Oct 11, 2024
Merged

[libc][math] Improve performance of double precision trig functions. #111793

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
fee360a
[libc][math] Improve performance of double precision trig functions.
lntue Oct 10, 2024
2f10f9c
Remove excessive `NO_FMA`.
lntue Oct 10, 2024
6c1cc4c
Fix formatting.
lntue Oct 10, 2024
cbced2b
Address comments.
lntue Oct 10, 2024
ed5f066
Remove extra ; in comments.
lntue Oct 10, 2024
56d5c65
Fix quick_mult.
lntue Oct 10, 2024
4450590
Fix exact_mult call in pow.cpp.
lntue Oct 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
43 changes: 26 additions & 17 deletions libc/src/__support/FPUtil/double_double.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@
namespace LIBC_NAMESPACE_DECL {
namespace fputil {

#define DEFAULT_DOUBLE_SPLIT 27

using DoubleDouble = LIBC_NAMESPACE::NumberPair<double>;

// The output of Dekker's FastTwoSum algorithm is correct, i.e.:
Expand Down Expand Up @@ -61,7 +63,8 @@ LIBC_INLINE constexpr DoubleDouble add(const DoubleDouble &a, double b) {
// Zimmermann, P., "Note on the Veltkamp/Dekker Algorithms with Directed
// Roundings," https://inria.hal.science/hal-04480440.
// Default splitting constant = 2^ceil(prec(double)/2) + 1 = 2^27 + 1.
template <size_t N = 27> LIBC_INLINE constexpr DoubleDouble split(double a) {
template <size_t N = DEFAULT_DOUBLE_SPLIT>
LIBC_INLINE constexpr DoubleDouble split(double a) {
DoubleDouble r{0.0, 0.0};
// CN = 2^N.
constexpr double CN = static_cast<double>(1 << N);
Expand All @@ -73,14 +76,30 @@ template <size_t N = 27> LIBC_INLINE constexpr DoubleDouble split(double a) {
return r;
}

// Helper for non-fma exact mult where the first number is already split.
template <size_t SPLIT_B = DEFAULT_DOUBLE_SPLIT>
LIBC_INLINE DoubleDouble exact_mult(const DoubleDouble &as, double a,
double b) {
DoubleDouble bs = split<SPLIT_B>(b);
DoubleDouble r{0.0, 0.0};

r.hi = a * b;
double t1 = as.hi * bs.hi - r.hi;
double t2 = as.hi * bs.lo + t1;
double t3 = as.lo * bs.hi + t2;
r.lo = as.lo * bs.lo + t3;

return r;
}

// Note: When FMA instruction is not available, the `exact_mult` function is
// only correct for round-to-nearest mode. See:
// Zimmermann, P., "Note on the Veltkamp/Dekker Algorithms with Directed
// Roundings," https://inria.hal.science/hal-04480440.
// Using Theorem 1 in the paper above, without FMA instruction, if we restrict
// the generated constants to precision <= 51, and splitting it by 2^28 + 1,
// then a * b = r.hi + r.lo is exact for all rounding modes.
template <bool NO_FMA_ALL_ROUNDINGS = false>
template <size_t SPLIT_B = 27>
Copy link
Member

@nickdesaulniers nickdesaulniers Oct 15, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

template <size_t SPLIT_B = DEFAULT_DOUBLE_SPLIT>

LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
DoubleDouble r{0.0, 0.0};

Expand All @@ -90,18 +109,8 @@ LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
#else
// Dekker's Product.
DoubleDouble as = split(a);
DoubleDouble bs;

if constexpr (NO_FMA_ALL_ROUNDINGS)
bs = split<28>(b);
else
bs = split(b);

r.hi = a * b;
double t1 = as.hi * bs.hi - r.hi;
double t2 = as.hi * bs.lo + t1;
double t3 = as.lo * bs.hi + t2;
r.lo = as.lo * bs.lo + t3;
r = exact_mult<SPLIT_B>(as, a, b);
#endif // LIBC_TARGET_CPU_HAS_FMA

return r;
Expand All @@ -113,10 +122,10 @@ LIBC_INLINE DoubleDouble quick_mult(double a, const DoubleDouble &b) {
return r;
}

template <bool NO_FMA_ALL_ROUNDINGS = false>
template <size_t SPLIT_B = 27>
Copy link
Member

@nickdesaulniers nickdesaulniers Oct 15, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

template <size_t SPLIT_B = DEFAULT_DOUBLE_SPLIT>

LIBC_INLINE DoubleDouble quick_mult(const DoubleDouble &a,
const DoubleDouble &b) {
DoubleDouble r = exact_mult<NO_FMA_ALL_ROUNDINGS>(a.hi, b.hi);
DoubleDouble r = exact_mult<SPLIT_B>(a.hi, b.hi);
double t1 = multiply_add(a.hi, b.lo, r.lo);
double t2 = multiply_add(a.lo, b.hi, t1);
r.lo = t2;
Expand Down Expand Up @@ -157,8 +166,8 @@ LIBC_INLINE DoubleDouble div(const DoubleDouble &a, const DoubleDouble &b) {
double e_hi = fputil::multiply_add(b.hi, -r.hi, a.hi);
double e_lo = fputil::multiply_add(b.lo, -r.hi, a.lo);
#else
DoubleDouble b_hi_r_hi = fputil::exact_mult</*NO_FMA=*/true>(b.hi, -r.hi);
DoubleDouble b_lo_r_hi = fputil::exact_mult</*NO_FMA=*/true>(b.lo, -r.hi);
DoubleDouble b_hi_r_hi = fputil::exact_mult(b.hi, -r.hi);
DoubleDouble b_lo_r_hi = fputil::exact_mult(b.lo, -r.hi);
double e_hi = (a.hi + b_hi_r_hi.hi) + b_hi_r_hi.lo;
double e_lo = (a.lo + b_lo_r_hi.hi) + b_lo_r_hi.lo;
#endif // LIBC_TARGET_CPU_HAS_FMA
Expand Down
10 changes: 10 additions & 0 deletions libc/src/__support/macros/optimization.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,16 @@ LIBC_INLINE constexpr bool expects_bool_condition(T value, T expected) {

#ifndef LIBC_MATH
#define LIBC_MATH 0
#else

#if (LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS)
#define LIBC_MATH_HAS_SKIP_ACCURATE_PASS
#endif

#if (LIBC_MATH & LIBC_MATH_SMALL_TABLES)
#define LIBC_MATH_HAS_SMALL_TABLES
#endif

#endif // LIBC_MATH

#endif // LLVM_LIBC_SRC___SUPPORT_MACROS_OPTIMIZATION_H
117 changes: 51 additions & 66 deletions libc/src/math/generic/cos.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,17 +17,14 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
#include "src/math/generic/range_reduction_double_common.h"
#include "src/math/generic/sincos_eval.h"

// TODO: We might be able to improve the performance of large range reduction of
// non-FMA targets further by operating directly on 25-bit chunks of 128/pi and
// pre-split SIN_K_PI_OVER_128, but that might double the memory footprint of
// those lookup table.
#include "range_reduction_double_common.h"

#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0)
#define LIBC_MATH_COS_SKIP_ACCURATE_PASS
#endif
#ifdef LIBC_TARGET_CPU_HAS_FMA
#include "range_reduction_double_fma.h"
#else
#include "range_reduction_double_nofma.h"
#endif // LIBC_TARGET_CPU_HAS_FMA

namespace LIBC_NAMESPACE_DECL {

Expand All @@ -42,22 +39,29 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) {

DoubleDouble y;
unsigned k;
generic::LargeRangeReduction<NO_FMA> range_reduction_large{};
LargeRangeReduction range_reduction_large{};

// |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA)
// |x| < 2^16.
if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) {
// |x| < 2^-27
if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) {
// Signed zeros.
if (LIBC_UNLIKELY(x == 0.0))
return 1.0;

// For |x| < 2^-27, |cos(x) - 1| < |x|^2/2 < 2^-54 = ulp(1 - 2^-53)/2.
return fputil::round_result_slightly_down(1.0);
// |x| < 2^-7
if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) {
// |x| < 2^-27
if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) {
// Signed zeros.
if (LIBC_UNLIKELY(x == 0.0))
return 1.0;

// For |x| < 2^-27, |cos(x) - 1| < |x|^2/2 < 2^-54 = ulp(1 - 2^-53)/2.
return fputil::round_result_slightly_down(1.0);
}
// No range reduction needed.
k = 0;
y.lo = 0.0;
y.hi = x;
} else {
// Small range reduction.
k = range_reduction_small(x, y);
}

// // Small range reduction.
k = range_reduction_small(x, y);
} else {
// Inf or NaN
if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) {
Expand All @@ -70,70 +74,51 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) {
}

// Large range reduction.
k = range_reduction_large.compute_high_part(x);
y = range_reduction_large.fast();
k = range_reduction_large.fast(x, y);
}

DoubleDouble sin_y, cos_y;

generic::sincos_eval(y, sin_y, cos_y);
[[maybe_unused]] double err = generic::sincos_eval(y, sin_y, cos_y);

// Look up sin(k * pi/128) and cos(k * pi/128)
// Memory saving versions:

// Use 128-entry table instead:
// DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 127];
// uint64_t sin_s = static_cast<uint64_t>((k + 128) & 128) << (63 - 7);
// sin_k.hi = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val();
// sin_k.lo = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val();
// DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 127];
// uint64_t cos_s = static_cast<uint64_t>((k + 64) & 128) << (63 - 7);
// cos_k.hi = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val();
// cos_k.lo = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val();

// Use 64-entry table instead:
// auto get_idx_dd = [](unsigned kk) -> DoubleDouble {
// unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
// DoubleDouble ans = SIN_K_PI_OVER_128[idx];
// if (kk & 128) {
// ans.hi = -ans.hi;
// ans.lo = -ans.lo;
// }
// return ans;
// };
// DoubleDouble sin_k = get_idx_dd(k + 128);
// DoubleDouble cos_k = get_idx_dd(k + 64);

#ifdef LIBC_MATH_HAS_SMALL_TABLES
// Memory saving versions. Use 65-entry table.
auto get_idx_dd = [](unsigned kk) -> DoubleDouble {
unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
DoubleDouble ans = SIN_K_PI_OVER_128[idx];
if (kk & 128) {
ans.hi = -ans.hi;
ans.lo = -ans.lo;
}
return ans;
};
DoubleDouble sin_k = get_idx_dd(k + 128);
DoubleDouble cos_k = get_idx_dd(k + 64);
#else
// Fast look up version, but needs 256-entry table.
// -sin(k * pi/128) = sin((k + 128) * pi/128)
// cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
DoubleDouble msin_k = SIN_K_PI_OVER_128[(k + 128) & 255];
DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255];
#endif // LIBC_MATH_HAS_SMALL_TABLES

// After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128).
// So k is an integer and -pi / 256 <= y <= pi / 256.
// Then cos(x) = cos((k * pi/128 + y)
// = cos(y) * cos(k*pi/128) - sin(y) * sin(k*pi/128)
DoubleDouble cos_k_cos_y = fputil::quick_mult<NO_FMA>(cos_y, cos_k);
DoubleDouble msin_k_sin_y = fputil::quick_mult<NO_FMA>(sin_y, msin_k);
DoubleDouble cos_k_cos_y = fputil::quick_mult(cos_y, cos_k);
DoubleDouble msin_k_sin_y = fputil::quick_mult(sin_y, msin_k);

DoubleDouble rr = fputil::exact_add<false>(cos_k_cos_y.hi, msin_k_sin_y.hi);
rr.lo += msin_k_sin_y.lo + cos_k_cos_y.lo;

#ifdef LIBC_MATH_COS_SKIP_ACCURATE_PASS
#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
return rr.hi + rr.lo;
#else

// Accurate test and pass for correctly rounded implementation.
#ifdef LIBC_TARGET_CPU_HAS_FMA
constexpr double ERR = 0x1.0p-70;
#else
// TODO: Improve non-FMA fast pass accuracy.
constexpr double ERR = 0x1.0p-66;
#endif // LIBC_TARGET_CPU_HAS_FMA

double rlp = rr.lo + ERR;
double rlm = rr.lo - ERR;
double rlp = rr.lo + err;
double rlm = rr.lo - err;

double r_upper = rr.hi + rlp; // (rr.lo + ERR);
double r_lower = rr.hi + rlm; // (rr.lo - ERR);
Expand All @@ -144,15 +129,15 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) {

Float128 u_f128, sin_u, cos_u;
if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT))
u_f128 = generic::range_reduction_small_f128(x);
u_f128 = range_reduction_small_f128(x);
else
u_f128 = range_reduction_large.accurate();

generic::sincos_eval(u_f128, sin_u, cos_u);

auto get_sin_k = [](unsigned kk) -> Float128 {
unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
Float128 ans = generic::SIN_K_PI_OVER_128_F128[idx];
Float128 ans = SIN_K_PI_OVER_128_F128[idx];
if (kk & 128)
ans.sign = Sign::NEG;
return ans;
Expand All @@ -172,7 +157,7 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) {
// https://github.com/llvm/llvm-project/issues/96452.

return static_cast<double>(r);
#endif // !LIBC_MATH_COS_SKIP_ACCURATE_PASS
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}

} // namespace LIBC_NAMESPACE_DECL
2 changes: 1 addition & 1 deletion libc/src/math/generic/pow.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -398,7 +398,7 @@ LLVM_LIBC_FUNCTION(double, pow, (double x, double y)) {
#else
double c = FPBits(m_x.uintval() & 0x3fff'e000'0000'0000).get_val();
dx = fputil::multiply_add(RD[idx_x], m_x.get_val() - c, CD[idx_x]); // Exact
dx_c0 = fputil::exact_mult<true>(COEFFS[0], dx);
dx_c0 = fputil::exact_mult<28>(dx, COEFFS[0]); // Exact
#endif // LIBC_TARGET_CPU_HAS_FMA

double dx2 = dx * dx;
Expand Down
Loading
Loading