Move primality test to Field_Zp::init

Throw exception when not prime
Add tests

7 files changed, 43 insertions, 26 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
diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx
index 28fbe3af..adc40499 100644
--- a/src/python/gudhi/cubical_complex.pyx
+++ b/src/python/gudhi/cubical_complex.pyx
@@ -35,7 +35,7 @@ cdef extern from "Cubical_complex_interface.h" namespace "Gudhi":
 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
     cdef cppclass Cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Cubical_complex::Cubical_complex_interface<>>":
         Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max) nogil
-        void compute_persistence(int homology_coeff_field, double min_persistence) nogil
+        void compute_persistence(int homology_coeff_field, double min_persistence) nogil except+
         vector[pair[int, pair[double, double]]] get_persistence() nogil
         vector[vector[int]] cofaces_of_cubical_persistence_pairs() nogil
         vector[int] betti_numbers() nogil
diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx
index d353d2af..0eaa5867 100644
--- a/src/python/gudhi/periodic_cubical_complex.pyx
+++ b/src/python/gudhi/periodic_cubical_complex.pyx
@@ -32,7 +32,7 @@ cdef extern from "Cubical_complex_interface.h" namespace "Gudhi":
 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
     cdef cppclass Periodic_cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Cubical_complex::Cubical_complex_interface<Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double>>>":
         Periodic_cubical_complex_persistence_interface(Periodic_cubical_complex_base_interface * st, bool persistence_dim_max) nogil
-        void compute_persistence(int homology_coeff_field, double min_persistence) nogil
+        void compute_persistence(int homology_coeff_field, double min_persistence) nogil except +
         vector[pair[int, pair[double, double]]] get_persistence() nogil
         vector[vector[int]] cofaces_of_cubical_persistence_pairs() nogil
         vector[int] betti_numbers() nogil
diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 3b8ea4f9..006a24ed 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -78,7 +78,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
     cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>":
         Simplex_tree_persistence_interface(Simplex_tree_interface_full_featured * st, bool persistence_dim_max) nogil
-        void compute_persistence(int homology_coeff_field, double min_persistence) nogil
+        void compute_persistence(int homology_coeff_field, double min_persistence) nogil except +
         vector[pair[int, pair[double, double]]] get_persistence() nogil
         vector[int] betti_numbers() nogil
         vector[int] persistent_betti_numbers(double from_value, double to_value) nogil
diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h
index 6877f190..e5a3dfba 100644
--- a/src/python/include/Persistent_cohomology_interface.h
+++ b/src/python/include/Persistent_cohomology_interface.h
@@ -43,21 +43,6 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol
     }
   };
 
-  bool is_prime(int n) {
-    // Primality test using 6k+-1 optimization.
-    if (n <= 3)
-      return n > 1;
-    if ((n % 2 == 0) || (n % 3 == 0))
-      return false;
-    int i = 5;
-    while (i*i <= n) {
-      if ((n % i == 0) || (n % (i + 2) == 0))
-        return false;
-      i += 6;
-    }
-    return true;
-  }
-
  public:
   Persistent_cohomology_interface(FilteredComplex* stptr, bool persistence_dim_max=false)
       : Base(*stptr, persistence_dim_max),
@@ -65,11 +50,6 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol
 
   // TODO: move to the constructors?
   void compute_persistence(int homology_coeff_field, double min_persistence) {
-    // Check that homology_coeff_field is a prime number (including not null)
-    if(!is_prime(homology_coeff_field)) {
-        std::cerr << "Warning: The persistence was not computed ; homology_coeff_field must be a prime number";
-        return;
-    }
     Base::init_coefficients(homology_coeff_field);
     Base::compute_persistent_cohomology(min_persistence);
   }