[libc-commits] [libc] [libc][math] Implement an integer-only version of double precision sin and cos with 1 ULP errors. (PR #184752)

[Simon Tatham via libc-commits]

================
@@ -0,0 +1,274 @@
+//===-- 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/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math_extras.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace integer_only {
+
+struct Frac128 : public UInt<128> {
+  using UInt<128>::UInt;
+
+  constexpr Frac128 operator~() const {
+    Frac128 r;
+    r.val[0] = ~val[0];
+    r.val[1] = ~val[1];
+    return r;
+  }
+
+  constexpr Frac128 operator+(const Frac128 &other) const {
+    UInt<128> r = UInt<128>(*this) + (UInt<128>(other));
+    return Frac128(r.val);
+  }
+
+  constexpr Frac128 operator-(const Frac128 &other) const {
+    UInt<128> r = UInt<128>(*this) - (UInt<128>(other));
+    return Frac128(r.val);
+  }
+
+  constexpr Frac128 operator*(const Frac128 &other) const {
+    UInt<128> r = UInt<128>::quick_mul_hi(UInt<128>(other));
+    return Frac128(r.val);
+  }
+};
+
+// 1280 + 64 bits of 2/pi, printed using MPFR.
+// Notice that if we store from the highest bytes to lowest bytes, it is
+// essentially having 2/pi in big-endian.  On the other hand, uint64_t type
+// that will be used for computations later are in little-endian.  So a few
----------------
statham-arm wrote:

Not necessarily! Big-endian architectures are still supported by LLVM. Arm can be BE, for example. I guess this comment has an implicit "... on the most usual architectures"?

As well as bit-reversal, another question is whether the architecture allows unaligned loads. If it doesn't, then even if the table endianness matches the CPU endianness, you still have to do multiple loads and bit-twiddling.

If performance of this part is a serious concern, perhaps it would be easiest to just store a sequence of uint64_t, and do two 64-bit loads followed by `(a << this) | (b >> that)`. An extra advantage is that then you get to index into 2/π with bit rather than byte precision.

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