Move primality test to Field_Zp::init

Throw exception when not prime
Add tests

3 files changed, 40 insertions, 3 deletions

diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt
index c68c6524..3e7e9369 100644
--- a/src/Persistent_cohomology/example/CMakeLists.txt
+++ b/src/Persistent_cohomology/example/CMakeLists.txt
@@ -11,7 +11,7 @@ if (TBB_FOUND)
   target_link_libraries(persistence_from_simple_simplex_tree ${TBB_LIBRARIES})
 endif()
 add_test(NAME Persistent_cohomology_example_from_simple_simplex_tree COMMAND $<TARGET_FILE:persistence_from_simple_simplex_tree>
-    "1" "0")
+    "2" "0")
 
 if(TARGET Boost::program_options)
   add_executable(rips_persistence_step_by_step rips_persistence_step_by_step.cpp)
diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h
index 0673625c..4bfd95c0 100644
--- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h
+++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h
@@ -13,6 +13,8 @@
 
 #include <utility>
 #include <vector>
+#include <stdexcept>
+#include <cmath>
 
 namespace Gudhi {
 
@@ -34,15 +36,30 @@ class Field_Zp {
 
   void init(int charac) {
     assert(charac > 0);  // division by zero + non negative values
+
     Prime = charac;
+
+    // Check for primality
+    if ((Prime == 0) || (Prime == 1) || ((Prime > 3) && ((Prime % 2 == 0) || (Prime % 3 == 0))))
+        throw std::invalid_argument("homology_coeff_field must be a prime number");
+
     inverse_.clear();
     inverse_.reserve(charac);
     inverse_.push_back(0);
     for (int i = 1; i < Prime; ++i) {
       int inv = 1;
-      while (((inv * i) % Prime) != 1)
+      int mult = inv * i;
+      while ( (mult % Prime) != 1) {
         ++inv;
+        if(mult == Prime)
+            throw std::invalid_argument("homology_coeff_field must be a prime number");
+        mult = inv * i;
+      }
       inverse_.push_back(inv);
+      if ( (i <= std::sqrt(Prime)) && (((i-5)%6) == 0) ) {
+        if ((Prime % i == 0) || (Prime % (i + 2) == 0))
+            throw std::invalid_argument("homology_coeff_field must be a prime number");
+      }
     }
   }
 
diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp
index fe3f8517..9559b842 100644
--- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp
+++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp
@@ -146,9 +146,14 @@ void test_rips_persistence_in_dimension(int dimension) {
   std::clog << "str_rips_persistence=" << str_rips_persistence << std::endl;
 }
 
+BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_0 )
+{
+  BOOST_CHECK_THROW(test_rips_persistence_in_dimension(0), std::invalid_argument);
+}
+
 BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_1 )
 {
-  test_rips_persistence_in_dimension(1);
+  BOOST_CHECK_THROW(test_rips_persistence_in_dimension(1), std::invalid_argument);
 }
 
 BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_2 )
@@ -161,11 +166,26 @@ BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_3 )
   test_rips_persistence_in_dimension(3);
 }
 
+BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_4 )
+{
+  BOOST_CHECK_THROW(test_rips_persistence_in_dimension(4), std::invalid_argument);
+}
+
 BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_5 )
 {
   test_rips_persistence_in_dimension(5);
 }
 
+BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_11 )
+{
+  test_rips_persistence_in_dimension(11);
+}
+
+BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_13 )
+{
+  test_rips_persistence_in_dimension(13);
+}
+
 // TODO(VR): not working from 6
 // std::string str_rips_persistence = test_rips_persistence(6, 0);
 // TODO(VR): division by zero