Merge pull request #348 from bluescarni/pr/vsop2013_fix

VSOP2013 fix

.clang-format

@@ -56,6 +56,6 @@ SpacesInContainerLiterals: false
 SpacesInCStyleCastParentheses: false
 SpacesInParentheses: false
 SpacesInSquareBrackets: false
-Standard:        Cpp11
+Standard:        Auto
 TabWidth:        4
 UseTab:          Never

.clang-tidy

@@ -1,5 +1,5 @@
 ---
-Checks:          'clang-diagnostic-*,clang-analyzer-*,*,-hicpp-uppercase-literal-suffix,-readability-uppercase-literal-suffix,-modernize-use-trailing-return-type,-readability-named-parameter,-hicpp-named-parameter,-cppcoreguidelines-pro-bounds-array-to-pointer-decay,-hicpp-no-array-decay,-llvm-header-guard,-cppcoreguidelines-macro-usage,-google-runtime-references,-readability-isolate-declaration,-fuchsia-default-arguments-calls,-fuchsia-overloaded-operator,-fuchsia-default-arguments-declarations,-readability-else-after-return,-google-runtime-int,-hicpp-signed-bitwise,-cert-dcl21-cpp,-cppcoreguidelines-avoid-magic-numbers,-readability-magic-numbers,-cppcoreguidelines-pro-bounds-pointer-arithmetic,-cppcoreguidelines-pro-type-reinterpret-cast,-cppcoreguidelines-avoid-c-arrays,-hicpp-avoid-c-arrays,-modernize-avoid-c-arrays,-modernize-use-transparent-functors,-cert-dcl16-c,-cppcoreguidelines-pro-type-union-access,-bugprone-branch-clone,-fuchsia-statically-constructed-objects,-cppcoreguidelines-pro-bounds-constant-array-index,-readability-static-accessed-through-instance,-cppcoreguidelines-pro-type-vararg,-hicpp-vararg,-llvmlibc-restrict-system-libc-headers,-llvmlibc-callee-namespace,-llvmlibc-implementation-in-namespace,-llvm-else-after-return,-fuchsia-trailing-return,-readability-identifier-length,-altera-unroll-loops,-altera-id-dependent-backward-branch,-altera-struct-pack-align,-performance-no-int-to-ptr,-readability-function-cognitive-complexity'
+Checks:          'clang-diagnostic-*,clang-analyzer-*,*,-hicpp-uppercase-literal-suffix,-readability-uppercase-literal-suffix,-modernize-use-trailing-return-type,-readability-named-parameter,-hicpp-named-parameter,-cppcoreguidelines-pro-bounds-array-to-pointer-decay,-hicpp-no-array-decay,-llvm-header-guard,-cppcoreguidelines-macro-usage,-google-runtime-references,-readability-isolate-declaration,-fuchsia-default-arguments-calls,-fuchsia-overloaded-operator,-fuchsia-default-arguments-declarations,-readability-else-after-return,-google-runtime-int,-hicpp-signed-bitwise,-cert-dcl21-cpp,-cppcoreguidelines-avoid-magic-numbers,-readability-magic-numbers,-cppcoreguidelines-pro-bounds-pointer-arithmetic,-cppcoreguidelines-pro-type-reinterpret-cast,-cppcoreguidelines-avoid-c-arrays,-hicpp-avoid-c-arrays,-modernize-avoid-c-arrays,-modernize-use-transparent-functors,-cert-dcl16-c,-cppcoreguidelines-pro-type-union-access,-bugprone-branch-clone,-fuchsia-statically-constructed-objects,-cppcoreguidelines-pro-bounds-constant-array-index,-readability-static-accessed-through-instance,-cppcoreguidelines-pro-type-vararg,-hicpp-vararg,-llvmlibc-restrict-system-libc-headers,-llvmlibc-callee-namespace,-llvmlibc-implementation-in-namespace,-llvm-else-after-return,-fuchsia-trailing-return,-readability-identifier-length,-altera-unroll-loops,-altera-id-dependent-backward-branch,-altera-struct-pack-align,-performance-no-int-to-ptr,-readability-function-cognitive-complexity,-misc-include-cleaner,-llvmlibc-inline-function-decl'
 WarningsAsErrors: '*'
 AnalyzeTemporaryDtors: false
 FormatStyle:     none

.github/workflows/gha_ci.yml

@@ -114,6 +114,12 @@ jobs:
           fetch-depth: 2
       - name: Build
         run: bash tools/gha_conda_coverage.sh
+  conda_llvm16_asan:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3
+      - name: Build
+        run: bash tools/gha_llvm16_conda_asan.sh
   conda_llvm15_asan:
     runs-on: ubuntu-latest
     steps:

CMakeLists.txt

@@ -11,7 +11,7 @@ if(NOT CMAKE_BUILD_TYPE)
 	FORCE)
 endif()
 
-project(heyoka VERSION 2.0.0 LANGUAGES CXX C)
+project(heyoka VERSION 2.0.1 LANGUAGES CXX C)
 
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" "${CMAKE_CURRENT_SOURCE_DIR}/cmake/yacma")
 

doc/changelog.rst

@@ -1,6 +1,16 @@
 Changelog
 =========
 
+2.0.1 (unreleased)
+------------------
+
+Fix
+~~~
+
+- Fix orbital elements singularity when using the VSOP2013
+  theory at low precision
+  (`#348 <https://github.com/bluescarni/heyoka/pull/348>`__).
+
 2.0.0 (2023-09-22)
 ------------------
 

include/heyoka/expression.hpp

@@ -20,11 +20,9 @@
 #include <optional>
 #include <ostream>
 #include <sstream>
-#include <stdexcept>
 #include <string>
 #include <type_traits>
 #include <unordered_map>
-#include <unordered_set>
 #include <utility>
 #include <variant>
 #include <vector>
@@ -79,7 +77,7 @@ class HEYOKA_DLL_PUBLIC expression
     template <typename Archive>
     void serialize(Archive &ar, unsigned)
     {
-        ar &m_value;
+        ar & m_value;
     }
 
 public:
@@ -138,9 +136,9 @@ HEYOKA_DLL_PUBLIC expression operator""_ldbl(unsigned long long);
 #if defined(HEYOKA_HAVE_REAL128)
 
 template <char... Chars>
-inline expression operator"" _f128()
+inline expression operator""_f128()
 {
-    return expression{mppp::literals::operator"" _rq<Chars...>()};
+    return expression{mppp::literals::operator""_rq < Chars... > ()};
 }
 
 #endif

src/llvm_state.cpp

@@ -364,7 +364,15 @@ struct llvm_state::jit {
         auto lljit = lljit_builder.create();
         // LCOV_EXCL_START
         if (!lljit) {
-            throw std::invalid_argument("Error creating an LLJIT object");
+            auto err = lljit.takeError();
+
+            std::string err_report;
+            llvm::raw_string_ostream ostr(err_report);
+
+            ostr << err;
+
+            throw std::invalid_argument(
+                fmt::format("Could not create an LLJIT object. The full error message is:\n{}", ostr.str()));
         }
         // LCOV_EXCL_STOP
         m_lljit = std::move(*lljit);
@@ -579,6 +587,7 @@ auto llvm_state_bc_to_module(const std::string &module_name, const std::string &
 
 } // namespace detail
 
+// NOLINTNEXTLINE(cppcoreguidelines-rvalue-reference-param-not-moved)
 llvm_state::llvm_state(std::tuple<std::string, unsigned, bool, bool, bool> &&tup)
     : m_jitter(std::make_unique<jit>()), m_opt_level(std::get<1>(tup)), m_fast_math(std::get<2>(tup)),
       m_force_avx512(std::get<3>(tup)), m_slp_vectorize(std::get<4>(tup)), m_module_name(std::move(std::get<0>(tup)))
@@ -1071,11 +1080,10 @@ void llvm_state::optimise()
     // NOTE: the reason for this inconsistency is that opt uses PB.parsePassPipeline()
     // (instead of PB.buildPerModuleDefaultPipeline()) to set up the optimisation
     // pipeline. Indeed, if we replace PB.buildPerModuleDefaultPipeline(ol) with
-    // PB.buildPerModuleDefaultPipeline(MPM, "default<O3>") (which
-    // corresponds to invoking "opt -passes='default<O3>'"), we do NOT need to set
-    // SLP vectorization on here to get the SLP vectorizer. Not sure if we should consider
-    // switching to this alternative way of setting up the optimisation pipeline
-    // in the future.
+    // PB.parsePassPipeline(MPM, "default<O3>") (which corresponds to invoking
+    // "opt -passes='default<O3>'"), we do NOT need to set SLP vectorization on
+    // here to get the SLP vectorizer. Not sure if we should consider switching to this
+    // alternative way of setting up the optimisation pipeline in the future.
     llvm::PipelineTuningOptions pto;
     pto.SLPVectorization = m_slp_vectorize;
     llvm::PassBuilder PB(m_jitter->m_tm.get(), pto);

src/model/vsop2013.cpp

@@ -324,6 +324,12 @@ std::vector<expression> vsop2013_cartesian_impl(std::uint32_t pl_idx, expression
                              qp = sqrt(qp_2);
                          });
 
+    // Check if qp is zero. This indicates a zero inclination orbit.
+    bool zero_inc = false;
+    if (const auto *num_ptr = std::get_if<number>(&qp.value())) {
+        zero_inc = is_zero(*num_ptr);
+    }
+
     // E, e, sqrt(1 - e**2), cos(i), sin(i), cos(Om), sin(Om), sin(E), cos(E)
     expression E, e, sqrt_1me2, ci, si, cOm, sOm, sin_E, cos_E;
     tbb::parallel_invoke(
@@ -337,12 +343,26 @@ std::vector<expression> vsop2013_cartesian_impl(std::uint32_t pl_idx, expression
             ci = 1_dbl - 2_dbl * qp_2;
             si = sqrt(1_dbl - ci * ci);
         },
-        [&]() { cOm = q / qp; }, [&]() { sOm = p / qp; });
+        // NOTE: for zero inclination, set conventionally Om to zero.
+        [&]() { cOm = zero_inc ? 1_dbl : q / qp; }, [&]() { sOm = zero_inc ? 0_dbl : p / qp; });
+
+    // Check for zero eccentricity.
+    bool zero_ecc = false;
+    // NOTE: not sure if this is possible at all, but better safe than sorry.
+    // Let's exclude it from code coverage checks for the time being.
+    // LCOV_EXCL_START
+    if (const auto *num_ptr = std::get_if<number>(&e.value())) {
+        zero_ecc = is_zero(*num_ptr);
+    }
+    // LCOV_EXCL_STOP
 
     // cos(om), sin(om), q1/a, q2/a.
     expression com, som, q1_a, q2_a;
-    tbb::parallel_invoke([&]() { com = (k * cOm + h * sOm) / e; }, [&]() { som = (h * cOm - k * sOm) / e; },
-                         [&]() { q1_a = cos_E - e; }, [&]() { q2_a = sqrt_1me2 * sin_E; });
+    tbb::parallel_invoke(
+        // NOTE: for zero eccentricity, set conventionally om to zero.
+        [&]() { com = zero_ecc ? 1_dbl : (k * cOm + h * sOm) / e; },
+        [&]() { som = zero_ecc ? 0_dbl : (h * cOm - k * sOm) / e; }, [&]() { q1_a = cos_E - e; },
+        [&]() { q2_a = sqrt_1me2 * sin_E; });
 
     // Prepare the entries of the rotation matrix, and a few auxiliary quantities.
     expression R00, R01, R10, R11, R20, R21, v_num, v_den;

src/taylor_02.cpp

@@ -975,6 +975,7 @@ std::pair<llvm::Value *, llvm::Type *> taylor_compute_jet_compact_mode(
         //
         // to pass the data necessary to the parallel workers.
         par_data_t = llvm::StructType::get(context, {builder.getInt32Ty(), ext_fp_ptr_t, ext_fp_ptr_t});
+        // NOLINTNEXTLINE(cppcoreguidelines-owning-memory)
         gl_par_data = new llvm::GlobalVariable(md, par_data_t, false, llvm::GlobalVariable::InternalLinkage,
                                                llvm::ConstantAggregateZero::get(par_data_t));
 

test/opt_checks.cpp

@@ -63,9 +63,10 @@ TEST_CASE("pow vect")
         REQUIRE(!boost::algorithm::contains(md_ir, "llvm.pow"));
     }
 
-#if LLVM_VERSION_MAJOR >= 16
+    // NOTE: no idea why, but in LLVM 17 it seems like there's no autovec
+    // happening at all in this specific case. This needs to be investigated.
+#if LLVM_VERSION_MAJOR == 16
 
-    // NOTE: LLVM16 is currently the version tested in the CI on arm64.
     if (tf.aarch64) {
         REQUIRE(!boost::algorithm::contains(md_ir, "llvm.pow"));
     }

test/vsop2013.cpp

@@ -14,6 +14,7 @@
 #include <vector>
 
 #include <heyoka/expression.hpp>
+#include <heyoka/llvm_state.hpp>
 #include <heyoka/model/vsop2013.hpp>
 #include <heyoka/taylor.hpp>
 
@@ -780,3 +781,26 @@ TEST_CASE("vsop2013 mus")
     REQUIRE(mus[3] == 8.9970116036316091182e-10);
     REQUIRE(mus[9] == 2.1886997654259696800e-12);
 }
+
+// Bug: low-precision solutions may have zero inclination,
+// which leads to singularities when converting to Cartesian.
+TEST_CASE("vsop2013 low prec zero inc")
+{
+    auto ex = vsop2013_cartesian(1, kw::time = "tm"_var, kw::thresh = 1e-1)[0];
+
+    llvm_state s;
+
+    add_cfunc<double>(s, "f", {ex});
+
+    s.compile();
+
+    auto *cf_ptr
+        = reinterpret_cast<void (*)(double *, const double *, const double *, const double *)>(s.jit_lookup("f"));
+
+    double out = 0, in = 0;
+
+    cf_ptr(&out, &in, nullptr, nullptr);
+
+    REQUIRE(!std::isnan(out));
+    REQUIRE(out != 0);
+}

tools/circleci_ubuntu_arm64.sh

@@ -14,7 +14,7 @@ wget https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforg
 export deps_dir=$HOME/local
 export PATH="$HOME/miniconda/bin:$PATH"
 bash miniconda.sh -b -p $HOME/miniconda
-mamba create -y -q -p $deps_dir cxx-compiler c-compiler cmake llvmdev tbb-devel tbb boost-cpp sleef xtensor xtensor-blas blas blas-devel fmt spdlog 'mppp>=0.27'
+mamba create -y -p $deps_dir cxx-compiler c-compiler cmake llvmdev tbb-devel tbb boost-cpp sleef xtensor xtensor-blas blas blas-devel fmt spdlog 'mppp>=0.27'
 source activate $deps_dir
 
 # Create the build dir and cd into it.

tools/gha_conda_release_static.sh

@@ -16,7 +16,7 @@ export PATH="$HOME/miniconda/bin:$PATH"
 bash miniconda.sh -b -p $HOME/miniconda
 conda config --add channels conda-forge
 conda config --set channel_priority strict
-conda create -y -q -p $deps_dir c-compiler cxx-compiler cmake llvmdev tbb-devel tbb boost-cpp 'mppp>=0.27' sleef xtensor xtensor-blas blas blas-devel fmt spdlog zlib
+conda create -y -q -p $deps_dir c-compiler cxx-compiler cmake 'llvmdev=16.*' tbb-devel tbb boost-cpp 'mppp>=0.27' sleef xtensor xtensor-blas blas blas-devel fmt spdlog zlib
 source activate $deps_dir
 
 # Create the build dir and cd into it.

tools/gha_llvm16_conda_asan.sh

@@ -0,0 +1,30 @@
+#!/usr/bin/env bash
+
+# Echo each command
+set -x
+
+# Exit on error.
+set -e
+
+# Core deps.
+sudo apt-get install wget
+
+# Install conda+deps.
+wget https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-Linux-x86_64.sh -O mambaforge.sh
+export deps_dir=$HOME/local
+export PATH="$HOME/mambaforge/bin:$PATH"
+bash mambaforge.sh -b -p $HOME/mambaforge
+mamba create -y -q -p $deps_dir c-compiler cxx-compiler cmake 'llvmdev=16.*' tbb-devel tbb boost-cpp 'mppp>=0.27' sleef xtensor xtensor-blas blas blas-devel 'fmt=9.*' spdlog
+source activate $deps_dir
+
+# Create the build dir and cd into it.
+mkdir build
+cd build
+
+# GCC build.
+cmake ../ -DCMAKE_PREFIX_PATH=$deps_dir -DCMAKE_BUILD_TYPE=Debug -DHEYOKA_BUILD_TESTS=yes -DHEYOKA_BUILD_TUTORIALS=ON -DHEYOKA_WITH_MPPP=yes -DHEYOKA_WITH_SLEEF=yes -DCMAKE_CXX_FLAGS="-fsanitize=address" -DBoost_NO_BOOST_CMAKE=ON
+make -j2 VERBOSE=1
+ctest -V -j2 -E vsop2013
+
+set +e
+set +x