Merge pull request #1039 from MClemot/ripser128

[RipsPersistenceDiagram] 128 bits indices option
topology-tool-kit · Jun 28, 2024 · 9a74f41 · 9a74f41
2 parents e34f787 + 9768876
commit 9a74f41
Show file tree

Hide file tree

Showing 3 changed files with 30 additions and 7 deletions.
diff --git a/core/base/ripsPersistenceDiagram/CMakeLists.txt b/core/base/ripsPersistenceDiagram/CMakeLists.txt
@@ -7,4 +7,11 @@ ttk_add_base_library(ripsPersistenceDiagram
     ripser.h
   DEPENDS
     common
-)
+)
+
+option(TTK_ENABLE_RIPSER_128BITS_IDS "Use 128 bits indices in Ripser" OFF)
+mark_as_advanced(TTK_ENABLE_RIPSER_128BITS_IDS)
+
+if (TTK_ENABLE_RIPSER_128BITS_IDS)
+    target_compile_definitions(ripsPersistenceDiagram PUBLIC TTK_ENABLE_RIPSER_128BITS_IDS)
+endif()
diff --git a/core/base/ripsPersistenceDiagram/ripser.h b/core/base/ripsPersistenceDiagram/ripser.h
@@ -18,7 +18,12 @@
 namespace ripser {
 
   using value_t = double;
+#if defined(TTK_ENABLE_RIPSER_128BITS_IDS) \
+  && (defined(__GNUC__) || defined(__clang__))
+  using index_t = __int128;
+#else
   using index_t = int64_t;
+#endif
   using coefficient_t = uint16_t;
 
   using simplex_t = std::vector<index_t>;

diff --git a/core/base/ripsPersistenceDiagram/ripserpy.cpp b/core/base/ripsPersistenceDiagram/ripserpy.cpp
@@ -91,15 +91,25 @@ static const std::string clear_line("\r\033[K");
 
 static const size_t num_coefficient_bits = 8;
 
-// 1L on windows is ALWAYS 32 bits, when on unix systems is pointer size
+#if defined(TTK_ENABLE_RIPSER_128BITS_IDS) \
+  && (defined(__GNUC__) || defined(__clang__))
 static const index_t max_simplex_index
   = (__int128(1) << (8 * sizeof(index_t) - 1 - num_coefficient_bits)) - 1;
 
-namespace std {
-  inline std::string to_string(__int128) {
-    return "";
-  }
-} // namespace std
+void check_overflow(index_t i) {
+  if
+#ifdef USE_COEFFICIENTS
+    (i > max_simplex_index)
+#else
+    (i < 0)
+#endif
+    throw std::overflow_error("simplex index " + std::to_string((uint64_t)i)
+                              + " in filtration is larger than maximum index");
+}
+#else
+// 1L on windows is ALWAYS 32 bits, when on unix systems is pointer size
+static const index_t max_simplex_index
+  = (uintptr_t(1) << (8 * sizeof(index_t) - 1 - num_coefficient_bits)) - 1;
 
 void check_overflow(index_t i) {
   if
@@ -112,6 +122,7 @@ void check_overflow(index_t i) {
                               + " in filtration is larger than maximum index "
                               + std::to_string(max_simplex_index));
 }
+#endif
 
 class binomial_coeff_table {
   /* Using flatten matrix */