[Support] Add mathematical constants

[evandro]

Summary

Add own version of the mathematical constants from the upcoming C++20 std::numbers.

[evandro]

Split the changes to the AMDGPU target in a separate patch.

[efriedma]

llvm/lib/CodeGen/SelectionDAG/SelectionDAGBuilder.cpp

5003

I don't like taking a 64-bit constant and truncating it to 32 bits. For most common mathematical constants, it probably produces the right result, but it's not correctly rounded in general.

llvm/lib/Transforms/Utils/SimplifyLibCalls.cpp

1943–1944

Do we actually have enough precision here if the type is long double?

[xbolva00]

Maybe take std::numbers implementation - llvm::numbers?

[evandro]

In D68257#1690565, @xbolva00 wrote:

Maybe take std::numbers implementation - llvm::numbers?

Please, elaborate.

llvm/lib/CodeGen/SelectionDAG/SelectionDAGBuilder.cpp

5003

Suggestions?

llvm/lib/Transforms/Utils/SimplifyLibCalls.cpp

1943–1944

The method is defined as static Constant *get(Type* Ty, double V). Likewise, the POSIX constant M_E is double precision only.

[efriedma]

llvm/lib/CodeGen/SelectionDAG/SelectionDAGBuilder.cpp

5003

You can have two constants: ln2 and ln2f, which are double and float respectively.

llvm/lib/Transforms/Utils/SimplifyLibCalls.cpp

1943–1944

I guess this isn't making the existing problem worse, but please leave a FIXME, at least.

[efriedma]

llvm/include/llvm/Support/MathExtras.h

69

Please use the correct number of digits (the smallest number of digits required to produce the correct double-precision result). Adding extra digits which don't actually affect the value is confusing.

[hiraditya]

LGTM as the comments have been addressed.

[efriedma]

llvm/include/llvm/Support/MathExtras.h

66

The correct value of sqrt(2) in double-precision is 1.4142135623730951.

And now I don't trust any of the other values...

71

Please mark the constants with "f", e.g. 2.718282f, so there isn't an extra rounding step.

llvm/lib/Transforms/Utils/SimplifyLibCalls.cpp

1944

I'm not sure this describes the issue correctly. You can specify a long double as a string, or raw bits. It can't interoperate with the native long double because that might have the wrong width.

[evandro]

llvm/include/llvm/Support/MathExtras.h

66

double has a precision of 15 or 16 significant digits. I don't understand why are you suggesting 17 significant digits when you asked to trim the precision down.

Besides, the reference I provided states that this value is 1.41421356237309505. Whether it's rounded to 1.4142135623730950 or 1.4142135623730951 is a bit moot, IMO.

llvm/lib/Transforms/Utils/SimplifyLibCalls.cpp

1944

What wording would you suggest, please?

[efriedma]

llvm/include/llvm/Support/MathExtras.h

66

I asked for "the smallest number of digits required to produce the correct double-precision result". This is what you get if, for example, you ask Python 2.7 or later to convert the value to a string with repr() (printf "import math\nprint(repr(math.sqrt(2)))" | python). 1.414213562373095 produces a value that's different by one ulp.

Yes, a one ulp difference is unlikely to matter for most uses, but if we're going to take the time to define these, we should define them correctly.

llvm/lib/Transforms/Utils/SimplifyLibCalls.cpp

1944

Maybe just FIXME: add more precise value of "e" for various long double types.

[evandro]

llvm/include/llvm/Support/MathExtras.h

66

You're assuming that Python is correct. bc says 1.41421356237309504880. glibc's math.h says 1.41421356237309504880 as well. And none of these is the same as your 1.4142135623730951.

As I said, the precision of double is 15 to 16 digits and of float, 6 to 7 digits. math.h defines them with 20 digits, which is probably an agreeable precision, yes? But I believe that we call all live with a difference of ±1ulp.

[evandro]

@efriedma, you had a point. I now trimmed the precision down to one digit short of seeing a change in the mantissa bits.

[evandro]

Added the hex FP constants for easier auditing.

[efriedma]

LGTM

llvm/include/llvm/Support/MathExtras.h

66

1.4142135623730951 is the shortest decimal representation that produces the same double-precision number as 0x1.6a09e667f3bcdP+0. It isn't "correct" in any other sense, sure.

A few extra digits is okay, I guess.

[evandro]

llvm/include/llvm/Support/MathExtras.h

66

Indeed, for some numbers more digits were necessary than for others. For uniformity's sake, I used the maximum from the number of digits required.

Thank you.

llvm/include/llvm/Support/MathExtras.h

 unsigned char _BitScanForward(unsigned long *_Index, unsigned long _Mask);
 unsigned char _BitScanForward64(unsigned long *_Index, unsigned __int64 _Mask);
 unsigned char _BitScanReverse(unsigned long *_Index, unsigned long _Mask);
 unsigned char _BitScanReverse64(unsigned long *_Index, unsigned __int64 _Mask);
 }
 #endif
 
 namespace llvm {
+
 /// The behavior an operation has on an input of 0.
 enum ZeroBehavior {
   /// The returned value is undefined.
   ZB_Undefined,
   /// The returned value is numeric_limits<T>::max()
   ZB_Max,
   /// The returned value is numeric_limits<T>::digits
   ZB_Width
 };
 
+/// Mathematical constants.
+namespace numbers {
+// TODO: Track C++20 std::numbers.
+// TODO: Favor using the hexadecimal FP constants (requires C++17).
+constexpr double e          = 2.7182818284590452354, // (0x1.5bf0a8b145749P+1) https://oeis.org/A001113
+                 egamma     = .57721566490153286061, // (0x1.2788cfc6fb619P-1) https://oeis.org/A001620
+                 ln2        = .69314718055994530942, // (0x1.62e42fefa39efP-1) https://oeis.org/A002162
+                 ln10       = 2.3025850929940456840, // (0x1.24bb1bbb55516P+1) https://oeis.org/A002392
+                 log2e      = 1.4426950408889634074, // (0x1.71547652b82feP+0)
+                 log10e     = .43429448190325182765, // (0x1.bcb7b1526e50eP-2)
+                 pi         = 3.1415926535897932385, // (0x1.921fb54442d18P+1) https://oeis.org/A000796
+                 inv_pi     = .31830988618379067154, // (0x1.45f306bc9c883P-2) https://oeis.org/A049541
+                 sqrtpi     = 1.7724538509055160273, // (0x1.c5bf891b4ef6bP+0) https://oeis.org/A002161
+                 inv_sqrtpi = .56418958354775628695, // (0x1.20dd750429b6dP-1) https://oeis.org/A087197
+                 sqrt2      = 1.4142135623730950488, // (0x1.6a09e667f3bcdP+0) https://oeis.org/A00219
+                 inv_sqrt2  = .70710678118654752440, // (0x1.6a09e667f3bcdP-1)
+                 sqrt3      = 1.7320508075688772935, // (0x1.bb67ae8584caaP+0) https://oeis.org/A002194
+                 inv_sqrt3  = .57735026918962576451, // (0x1.279a74590331cP-1)
+                 phi        = 1.6180339887498948482; // (0x1.9e3779b97f4a8P+0) https://oeis.org/A001622
+constexpr float ef          = 2.71828183F, // (0x1.5bf0a8P+1) https://oeis.org/A001113
+                egammaf     = .577215665F, // (0x1.2788d0P-1) https://oeis.org/A001620
+                ln2f        = .693147181F, // (0x1.62e430P-1) https://oeis.org/A002162
+                ln10f       = 2.30258509F, // (0x1.26bb1cP+1) https://oeis.org/A002392
+                log2ef      = 1.44269504F, // (0x1.715476P+0)
+                log10ef     = .434294482F, // (0x1.bcb7b2P-2)
+                pif         = 3.14159265F, // (0x1.921fb6P+1) https://oeis.org/A000796
+                inv_pif     = .318309886F, // (0x1.45f306P-2) https://oeis.org/A049541
+                sqrtpif     = 1.77245385F, // (0x1.c5bf8aP+0) https://oeis.org/A002161
+                inv_sqrtpif = .564189584F, // (0x1.20dd76P-1) https://oeis.org/A087197
+                sqrt2f      = 1.41421356F, // (0x1.6a09e6P+0) https://oeis.org/A002193
+                inv_sqrt2f  = .707106781F, // (0x1.6a09e6P-1)
+                sqrt3f      = 1.73205081F, // (0x1.bb67aeP+0) https://oeis.org/A002194
+                inv_sqrt3f  = .577350269F, // (0x1.279a74P-1)
+                phif        = 1.61803399F; // (0x1.9e377aP+0) https://oeis.org/A001622
+} // namespace numbers
+
 namespace detail {
 template <typename T, std::size_t SizeOfT> struct TrailingZerosCounter {
   static unsigned count(T Val, ZeroBehavior) {
     if (!Val)
       return std::numeric_limits<T>::digits;
     if (Val & 0x1)
       return 0;
 

llvm/lib/CodeGen/SelectionDAG/SelectionDAGBuilder.cpp

 static SDValue expandExp(const SDLoc &dl, SDValue Op, SelectionDAG &DAG,
                          const TargetLowering &TLI) {
   if (Op.getValueType() == MVT::f32 &&
       LimitFloatPrecision > 0 && LimitFloatPrecision <= 18) {
 
     // Put the exponent in the right bit position for later addition to the
     // final result:
     //
-    //   #define LOG2OFe 1.4426950f
-    //   t0 = Op * LOG2OFe
+    // t0 = Op * log2(e)
 
     // TODO: What fast-math-flags should be set here?
     SDValue t0 = DAG.getNode(ISD::FMUL, dl, MVT::f32, Op,
-                             getF32Constant(DAG, 0x3fb8aa3b, dl));
+                             DAG.getConstantFP(numbers::log2ef, dl, MVT::f32));
     return getLimitedPrecisionExp2(t0, dl, DAG);
   }
 
   // No special expansion.
   return DAG.getNode(ISD::FEXP, dl, Op.getValueType(), Op);
 }
 
 /// expandLog - Lower a log intrinsic. Handles the special sequences for
 /// limited-precision mode.
 static SDValue expandLog(const SDLoc &dl, SDValue Op, SelectionDAG &DAG,
                          const TargetLowering &TLI) {
   // TODO: What fast-math-flags should be set on the floating-point nodes?
 
   if (Op.getValueType() == MVT::f32 &&
       LimitFloatPrecision > 0 && LimitFloatPrecision <= 18) {
     SDValue Op1 = DAG.getNode(ISD::BITCAST, dl, MVT::i32, Op);
 
-    // Scale the exponent by log(2) [0.69314718f].
+    // Scale the exponent by log(2).
     SDValue Exp = GetExponent(DAG, Op1, TLI, dl);
-    SDValue LogOfExponent = DAG.getNode(ISD::FMUL, dl, MVT::f32, Exp,
-                                        getF32Constant(DAG, 0x3f317218, dl));
+    SDValue LogOfExponent =
+        DAG.getNode(ISD::FMUL, dl, MVT::f32, Exp,
+                    DAG.getConstantFP(numbers::ln2f, dl, MVT::f32));
 
     // Get the significand and build it into a floating-point number with
     // exponent of 1.
     SDValue X = GetSignificand(DAG, Op1, dl);
 
     SDValue LogOfMantissa;
     if (LimitFloatPrecision <= 6) {
       // For floating-point precision of 6:

llvm/lib/Transforms/Utils/SimplifyLibCalls.cpp

 using namespace PatternMatch;
 
 static cl::opt<bool>
     EnableUnsafeFPShrink("enable-double-float-shrink", cl::Hidden,
                          cl::init(false),
                          cl::desc("Enable unsafe double to float "
                                   "shrinking for math lib calls"));
 
-
 //===----------------------------------------------------------------------===//
 // Helper Functions
 //===----------------------------------------------------------------------===//
 
 static bool ignoreCallingConv(LibFunc Func) {
   return Func == LibFunc_abs || Func == LibFunc_labs ||
          Func == LibFunc_llabs || Func == LibFunc_strlen;
 }
if (ArgLb == PowLb || ArgID == Intrinsic::pow) {
     return MulY;
   }
   // log(exp{,2,10}(y)) -> y*log({e,2,10})
   // TODO: There is no exp10() intrinsic yet.
   else if (ArgLb == ExpLb || ArgLb == Exp2Lb || ArgLb == Exp10Lb ||
            ArgID == Intrinsic::exp || ArgID == Intrinsic::exp2) {
     Constant *Eul;
     if (ArgLb == ExpLb || ArgID == Intrinsic::exp)
-      // FIXME: The Euler number should be M_E, but it's place of definition
-      // is not quite standard.
-      Eul = ConstantFP::get(Log->getType(), 2.7182818284590452354);
+      // FIXME: Add more precise value of e for long double.
+      Eul = ConstantFP::get(Log->getType(), numbers::e);
     else if (ArgLb == Exp2Lb || ArgID == Intrinsic::exp2)
       Eul = ConstantFP::get(Log->getType(), 2.0);
     else
       Eul = ConstantFP::get(Log->getType(), 10.0);
     Value *LogE = Log->doesNotAccessMemory()
                       ? B.CreateCall(Intrinsic::getDeclaration(Mod, LogID, Ty),
                                      Eul, "log")
                       : emitUnaryFloatFnCall(Eul, LogNm, B, Attrs);