[libc-commits] [libc] [libc][math] Implement an integer-only version of double precision sin and cos with 1 ULP errors. (PR #184752)
Muhammad Bassiouni via libc-commits
libc-commits at lists.llvm.org
Mon Mar 9 22:46:50 PDT 2026
================
@@ -0,0 +1,271 @@
+//===-- Trig range reduction and evaluation using integer-only --*- 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___SUPPORT_MATH_SINCOS_INTEGER_UTILS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_SINCOS_INTEGER_UTILS_H
+
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/big_int.h"
+#include "src/__support/frac128.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math_extras.h"
+
+#undef LIBC_TARGET_IS_BIG_ENDIAN
+#if !defined(__BYTE_ORDER__) || !defined(__ORDER_LITTLE_ENDIAN__) || \
+ !defined(__ORDER_BIG_ENDIAN__)
+#define LIBC_TARGET_IS_BIG_ENDIAN 0
+#else
+#define LIBC_TARGET_IS_BIG_ENDIAN (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
+#endif // /LIBC_TARGET_IS_BIG_ENDIAN
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace integer_only {
+
+// 1280 + 64 bits of 2/pi, printed using MPFR.
+// We also add 8 more bytes to extend to all non-negative exponents.
+LIBC_INLINE_VAR constexpr unsigned TWO_OVER_PI_LENGTH = 1280 / 8 + 7;
+
+#if LIBC_TARGET_IS_BIG_ENDIAN
+LIBC_INLINE_VAR constexpr uint8_t TWO_OVER_PI[TWO_OVER_PI_LENGTH] = {
+ 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xA2, 0xF9, 0x83, 0x6E, 0x4E,
+ 0x44, 0x15, 0x29, 0xFC, 0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB,
+ 0x62, 0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63, 0xAB, 0xDE,
+ 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A, 0x42, 0x4D, 0xD2, 0xE0, 0x06,
+ 0x49, 0x2E, 0xEA, 0x09, 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1,
+ 0x29, 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44, 0x84, 0xE9,
+ 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41, 0x39, 0x91, 0xD6, 0x39, 0x83,
+ 0x53, 0x39, 0xF4, 0x9C, 0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F,
+ 0xF8, 0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11, 0x8B, 0x5A,
+ 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF, 0x27, 0xCB, 0x09, 0xB7, 0x4F,
+ 0x46, 0x3F, 0x66, 0x9E, 0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB,
+ 0xE5, 0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92, 0xEA, 0x6B,
+ 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08, 0x56, 0x03, 0x30, 0x46, 0xFC,
+ 0x7B, 0x6B, 0xAB, 0xF0, 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D,
+};
+#else // !LIBC_TARGET_IS_BIG_ENDIAN
+LIBC_INLINE_VAR constexpr uint8_t TWO_OVER_PI[TWO_OVER_PI_LENGTH] = {
+ 0x1D, 0x36, 0xF4, 0x9A, 0x20, 0xBC, 0xCF, 0xF0, 0xAB, 0x6B, 0x7B, 0xFC,
+ 0x46, 0x30, 0x03, 0x56, 0x08, 0x5D, 0x8D, 0x1F, 0xB1, 0x5F, 0xFB, 0x6B,
+ 0xEA, 0x92, 0x52, 0x8A, 0xF7, 0x39, 0x07, 0x3D, 0x7B, 0xF1, 0xE5, 0xEB,
+ 0xC7, 0xBA, 0x27, 0x75, 0x2D, 0xEA, 0x5F, 0x9E, 0x66, 0x3F, 0x46, 0x4F,
+ 0xB7, 0x09, 0xCB, 0x27, 0xCF, 0x7E, 0x36, 0x6D, 0x1F, 0x6D, 0x0A, 0x5A,
+ 0x8B, 0x11, 0x2F, 0xEF, 0x0F, 0x98, 0x05, 0xDE, 0xFF, 0x97, 0xF8, 0x1F,
+ 0x3B, 0x28, 0xF9, 0xBD, 0x8B, 0x5F, 0x84, 0x9C, 0xF4, 0x39, 0x53, 0x83,
+ 0x39, 0xD6, 0x91, 0x39, 0x41, 0x7E, 0x5F, 0xB4, 0x26, 0x70, 0x9C, 0xE9,
+ 0x84, 0x44, 0xBB, 0x2E, 0xF5, 0x35, 0x82, 0xE8, 0x3E, 0xA7, 0x29, 0xB1,
+ 0x1C, 0xEB, 0x1D, 0xFE, 0x1C, 0x92, 0xD1, 0x09, 0xEA, 0x2E, 0x49, 0x06,
+ 0xE0, 0xD2, 0x4D, 0x42, 0x3A, 0x6E, 0x24, 0xB7, 0x61, 0xC5, 0xBB, 0xDE,
+ 0xAB, 0x63, 0x51, 0xFE, 0x41, 0x90, 0x43, 0x3C, 0x99, 0x95, 0x62, 0xDB,
+ 0xC0, 0xDD, 0x34, 0xF5, 0xD1, 0x57, 0x27, 0xFC, 0x29, 0x15, 0x44, 0x4E,
+ 0x6E, 0x83, 0xF9, 0xA2, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+};
+#endif // LIBC_TARGET_IS_BIG_ENDIAN
+
+LIBC_INLINE_VAR constexpr Frac128 PI_OVER_2_M1({0x898c'c517'01b8'39a2,
+ 0x921f'b544'42d1'8469});
+
+// Perform range reduction mod pi/2
+//
+// Inputs:
+// x_u: explicit mantissa
+// x_e: biased exponent
+// Output:
+// k : round(x * 2/pi) mod 4
+// x_frac: |x - k * pi/2|
+// Return:
+// x_frac_is_neg.
+LIBC_INLINE constexpr bool trig_range_reduction(uint64_t x_u, unsigned x_e,
+ unsigned &k, Frac128 &x_frac) {
+ using FPBits = typename fputil::FPBits<double>;
+ bool x_frac_is_neg = false;
+ // We do multiplication x * (2/pi)
+ // Let T[i] be the i'th byte of 2/pi expansion:
+ // Then 2/pi = T[0] * 2^-8 + T[1] * 2^-16 + ...
+ // = sum_i T[i] * 2^(-8(i + 1))
+ // To be able to drop all T[j] * 2^(-8(j + 1)) for small j < i, we will want
+ // ulp(x) * lsb(T[i - 1] * 2^(-8 * i)) >= 4 = 2^2 (since 4 * pi/2 = 2*pi)
+ // So:
+ // 2^(e - 52) * 2^(-8 * i) >= 2^2
+ // Or equivalently,
+ // e - 54 - 8*i >= 0.
+ // Define:
+ // i = floor( (e - 54)/8 ),
+ // and let
+ // s = e - 54 - 8i >= 0.
+ // Since we store the mantissa of x, which is 53 bits long in a 64 bit
+ // integer, we have some wiggle room to shuffle the lsb of x.
+ // By shifting mantissa of x_u to the left by s, the lsb of x_u will be:
+ // 2^(e - 52 - s), for which, the product of lsb's is now exactly 4
+ // lsb(x_u) * 2^(-8 * i)) = 4.
+ // This will allow us to compute the full product:
+ // x_u * (T[i] * 2^(-8(i + 1)) + ... ) in exact fixed point.
+ // From the formula of i, in order for i >= 0, e >= 54. To support all the
+ // exponents e >= 0, we could add ceil(54 / 8) = 7 0x00 bytes and shift the
+ // index by 7.
+ unsigned e_num = x_e - FPBits::EXP_BIAS + 2; // e - 54 + 7*8
+ // With
+ // i = floor( (e - 54) / 8 ),
+ // the shifted-by-7 index is:
+ // j = i + 7 = floor( (e - 54) / 8 ) + 7
+ // Since the 64-bit integer chunk will be form by T[j] ... T[j + 7],
+ // and we store the table in the little-endian form, we will index to the
+ // lowest part of the 64-bit integer chunk, which is:
+ // idx = the index of the T[j + 7] part.
+ unsigned j = e_num >> 3;
+ unsigned idx =
+ LIBC_TARGET_IS_BIG_ENDIAN ? j : (TWO_OVER_PI_LENGTH - 1 - j - 7);
+ unsigned shift = e_num & 7; // s = e - 54 - 8*i
+ x_u <<= shift; // lsb(x_u) = 2^(e - 52 - s)
+ UInt<64> x_u64(x_u);
+ // Gather parts
+#if LIBC_TARGET_IS_BIG_ENDIAN
+ auto get_uint64 = [](const uint8_t *ptr) -> uint64_t {
+ return ptr[7] | (uint64_t(ptr[6]) << 8) | (uint64_t(ptr[5]) << 16) |
+ (uint64_t(ptr[4]) << 24) | (uint64_t(ptr[3]) << 32) |
+ (uint64_t(ptr[2]) << 40) | (uint64_t(ptr[1]) << 48) |
+ (uint64_t(ptr[0]) << 56);
+ };
+#else // !LIBC_TARGET_IS_BIG_ENDIAN
+ auto get_uint64 = [](const uint8_t *ptr) -> uint64_t {
+ return ptr[0] | (uint64_t(ptr[1]) << 8) | (uint64_t(ptr[2]) << 16) |
+ (uint64_t(ptr[3]) << 24) | (uint64_t(ptr[4]) << 32) |
+ (uint64_t(ptr[5]) << 40) | (uint64_t(ptr[6]) << 48) |
+ (uint64_t(ptr[7]) << 56);
+ };
+#endif // LIBC_TARGET_IS_BIG_ENDIAN
+ // lsb(c0) = 2^(-8i - 64)
+ uint64_t c0 = get_uint64(&TWO_OVER_PI[idx]);
+ // lsb(p0) = lsb(x_u) * lsb(c0)
+ // = 2^(e - 52 - s) * 2^(-8i - 64)
+ // = 2^(-62)
+ // msb(p0) = 2^(-62 + 63) = 2^1.
+ uint64_t p0 = x_u * c0;
+ // lsb(c1) = lsb(c0) * 2^-64 = 2^(-8i - 128)
+ // lsb(c2) = lsb(c1) * 2^-64 = 2^(-8i - 192)
+#if LIBC_TARGET_IS_BIG_ENDIAN
+ UInt<64> c1(get_uint64(&TWO_OVER_PI[idx + 8]));
+ UInt<64> c2(get_uint64(&TWO_OVER_PI[idx + 16]));
+#else // !LIBC_TARGET_IS_BIG_ENDIAN
+ UInt<64> c1(get_uint64(&TWO_OVER_PI[idx - 8]));
+ UInt<64> c2(get_uint64(&TWO_OVER_PI[idx - 16]));
+#endif // LIBC_TARGET_IS_BIG_ENDIAN
+ // lsb(p1) = lsb(x_u) * lsb(c1) = 2^(-62 - 64) = 2^-126
+ UInt<128> p1 = x_u64.ful_mul(c1);
+ // lsb(p2) = lsb(x_u) * lsb(c2) * 2^64 = 2^-126
+ UInt<128> p2(x_u64.quick_mul_hi(c2));
+ UInt<128> sum = p1 + p2;
+ sum.val[1] += p0;
+ // Get the highest 2 bits.
+ k = static_cast<unsigned>(sum.val[1] >> 62);
+ bool round_bit = sum.val[1] & 0x2000'0000'0000'0000;
+ // Shift so that the leading bit is 0.5.
+ sum <<= 2;
+ x_frac = Frac128(sum.val);
+ // Round to nearest k.
+ if (round_bit) {
+ // Flip the sign.
+ x_frac_is_neg = true;
+ ++k;
+ // Fast approximation of `1 - x_frac` with error = -lsb(x_frac) = -2^-128.
+ // Since in 2-complement, -x = ~x + lsb(x).
+ x_frac = ~x_frac;
+ }
+
+ // Perform multiplication x_frac * pi/2
+ x_frac = fputil::multiply_add(x_frac, PI_OVER_2_M1, x_frac);
+
+ return x_frac_is_neg;
+}
+
+// 128-bit fixed-point minimax polynomial approximation of sin(x) generated by
+// Sollya with:
+// > P = fpminimax(sin(x), [|1, 3, 5, 7, 9, 11, 13|], [|1, 128...|],
+// [0, pi/4], fixed);
+// > dirtyinfnorm( (sin(x) - P(x))/sin(x), [0, pi/4]);
+// 0x1.17a4...p-58
+// Storing absolute values of the coefficients.
+LIBC_INLINE_VAR constexpr Frac128 SIN_COEFF[] = {
+ Frac128({0x91b3'96a3'd5c5'fd6a, 0x2aaa'aaaa'aaaa'8ff2}), // x^3
+ Frac128({0x321f'bc0b'b8ca'f059, 0x0222'2222'221e'eac3}), // x^5
+ Frac128({0x36aa'355c'3311'996d, 0x000d'00d0'0cdf'8c9b}), // x^7
+ Frac128({0x0556'929e'ad60'7cb2, 0x0000'2e3b'c6ab'd75e}), // x^9
+ Frac128({0xa260'c74f'239d'd891, 0x0000'006b'9795'15a2}), // x^11
+ Frac128({0x4c97'758e'92ac'214c, 0x0000'0000'aec7'1a39}), // x^13
+};
+// 128-bit fixed-point minimax polynomial approximation of cos(x) generated by
+// Sollya with:
+// > P = fpminimax(cos(x), [|0, 2, 4, 6, 8, 10, 12|], [|1, 128...|],
+// [0, pi/4], fixed);
+// > dirtyinfnorm( (cos(x) - P(x))/cos(x), [0, pi/4]);
+// 0x1.269f...p-54
+// Storing absolute values of the coefficients.
+LIBC_INLINE_VAR constexpr Frac128 COS_COEFF[] = {
+ Frac128({0x56f6'2e74'b16e'5555, 0x7fff'ffff'fffe'4bfe}), // x^2
+ Frac128({0x860a'3e6c'cc50'e0d8, 0x0aaa'aaaa'aa77'5c33}), // x^4
+ Frac128({0xa87a'8f81'7440'7dd6, 0x005b'05b0'58fc'6fed}), // x^6
+ Frac128({0x84b2'76a3'c971'e7b8, 0x0001'a019'f80a'8ad5}), // x^8
+ Frac128({0x0082'310d'4e65'6b1f, 0x0000'049f'7cff'73d2}), // x^10
+ Frac128({0xed56'891e'f750'c7a9, 0x0000'0008'dc50'133d}), // x^12
+};
+
+// Compute sin(x) with relative errors ~ 2^-54.
+LIBC_INLINE constexpr double sin_eval(const Frac128 &x_frac, unsigned k,
----------------
bassiounix wrote:
I think it's not done!
https://github.com/llvm/llvm-project/pull/184752
More information about the libc-commits
mailing list