Merge pull request #452 from bluescarni/pr/new_llvm_math

New LLVM math primitives

include/heyoka/detail/llvm_helpers.hpp

@@ -157,6 +157,10 @@ HEYOKA_DLL_PUBLIC llvm::Value *llvm_ui_to_fp(llvm_state &, llvm::Value *, llvm::
 
 HEYOKA_DLL_PUBLIC llvm::Value *llvm_abs(llvm_state &, llvm::Value *);
 HEYOKA_DLL_PUBLIC llvm::Value *llvm_floor(llvm_state &, llvm::Value *);
+HEYOKA_DLL_PUBLIC llvm::Value *llvm_trunc(llvm_state &, llvm::Value *);
+
+HEYOKA_DLL_PUBLIC llvm::Value *llvm_is_finite(llvm_state &, llvm::Value *);
+HEYOKA_DLL_PUBLIC llvm::Value *llvm_is_natural(llvm_state &, llvm::Value *);
 
 HEYOKA_DLL_PUBLIC std::pair<llvm::Value *, llvm::Value *> llvm_sincos(llvm_state &, llvm::Value *);
 HEYOKA_DLL_PUBLIC llvm::Value *llvm_atan2(llvm_state &, llvm::Value *, llvm::Value *);

include/heyoka/detail/real_helpers.hpp

@@ -44,6 +44,7 @@ llvm::Value *llvm_real_fcmp_one(llvm_state &, llvm::Value *, llvm::Value *);
 llvm::Value *llvm_real_fnz(llvm_state &, llvm::Value *);
 llvm::Value *llvm_real_ui_to_fp(llvm_state &, llvm::Value *, llvm::Type *);
 llvm::Value *llvm_real_sgn(llvm_state &, llvm::Value *);
+llvm::Value *llvm_real_isfinite(llvm_state &, llvm::Value *);
 
 HEYOKA_DLL_PUBLIC mppp::real eps_from_prec(mpfr_prec_t);
 

src/detail/llvm_helpers.cpp

@@ -2347,6 +2347,69 @@ llvm::Value *llvm_floor(llvm_state &s, llvm::Value *x)
                           x);
 }
 
+// Trunc.
+llvm::Value *llvm_trunc(llvm_state &s, llvm::Value *x)
+{
+    return llvm_math_intr(s, "llvm.trunc",
+#if defined(HEYOKA_HAVE_REAL128)
+                          "truncq",
+#endif
+#if defined(HEYOKA_HAVE_REAL)
+                          "mpfr_trunc",
+#endif
+                          x);
+}
+
+// is_finite().
+llvm::Value *llvm_is_finite(llvm_state &s, llvm::Value *x)
+{
+    assert(x != nullptr);
+
+    auto *x_t = x->getType();
+
+    if (x_t->getScalarType()->isFloatingPointTy()) {
+        // Codegen +- inf.
+        auto *pinf = llvm_codegen(s, x_t, number{std::numeric_limits<double>::infinity()});
+        auto *minf = llvm_codegen(s, x_t, number{-std::numeric_limits<double>::infinity()});
+
+        // Check that if x is not +- inf or NaN.
+        auto *x_not_pinf = llvm_fcmp_one(s, x, pinf);
+        auto *x_not_minf = llvm_fcmp_one(s, x, minf);
+
+        // Put the conditions together and return.
+        return s.builder().CreateLogicalAnd(x_not_pinf, x_not_minf);
+#if defined(HEYOKA_HAVE_REAL)
+    } else if (llvm_is_real(x_t) != 0) {
+        return llvm_real_isfinite(s, x);
+#endif
+        // LCOV_EXCL_START
+    } else [[unlikely]] {
+        throw std::invalid_argument(fmt::format(
+            "Invalid type '{}' encountered in the LLVM implementation of is_finite()", llvm_type_name(x_t)));
+    }
+    // LCOV_EXCL_STOP
+}
+
+// Check if the input floating-point value is a natural number.
+llvm::Value *llvm_is_natural(llvm_state &s, llvm::Value *x)
+{
+    // Is x finite?
+    auto *x_finite = llvm_is_finite(s, x);
+
+    // Is x>=0?
+    auto *x_ge_0 = llvm_fcmp_oge(s, x, llvm_codegen(s, x->getType(), number{0.}));
+
+    // Is x an integral value?
+    auto *x_int = llvm_fcmp_oeq(s, x, llvm_trunc(s, x));
+
+    // Put the conditions together and return.
+    auto &bld = s.builder();
+    auto *ret = bld.CreateLogicalAnd(x_finite, x_ge_0);
+    ret = bld.CreateLogicalAnd(ret, x_int);
+
+    return ret;
+}
+
 // Add a function to count the number of sign changes in the coefficients
 // of a polynomial of degree n. The coefficients are SIMD vectors of size batch_size
 // and scalar type scal_t.

src/detail/real_helpers.cpp

@@ -611,6 +611,21 @@ llvm::Value *llvm_real_fnz(llvm_state &s, llvm::Value *x)
     return builder.CreateICmpEQ(ret, llvm::ConstantInt::getNullValue(ret->getType()));
 }
 
+llvm::Value *llvm_real_isfinite(llvm_state &s, llvm::Value *x)
+{
+    // LCOV_EXCL_START
+    assert(x != nullptr);
+    // LCOV_EXCL_STOP
+
+    auto &bld = s.builder();
+
+    // Check if x is an ordinary number.
+    auto *f = real_nary_cmp(s, x->getType(), "mpfr_number_p", 1u);
+    auto *ret = bld.CreateCall(f, x);
+
+    return ret;
+}
+
 // Convert the input unsigned integral value n to the real type fp_t.
 llvm::Value *llvm_real_ui_to_fp(llvm_state &s, llvm::Value *n, llvm::Type *fp_t)
 {

src/detail/vector_math.cpp

@@ -126,7 +126,7 @@ auto make_vf_map()
     // by sleef, on the assumption that usually sqrt() is implemented directly in hardware
     // and thus there's no need to go through sleef. This is certainly true for x86,
     // but I am not 100% sure for the other archs. Let's keep this in mind.
-    // NOTE: the same holds for things like abs() and floor().
+    // NOTE: the same holds for things like abs(), floor(), trunc(), etc.
 
     // Single-precision.
     add_vfinfo_sleef(retval, "llvm.sin.f32", "sin", "f");