Verification of perioidcity congruences is finished

Verification of periodicity congruences is finished and works ok.

.gitignore

@@ -1,5 +1,6 @@
 #*
 *~
+.*
 *.o
 GPATH
 GRTAGS

Makefile

@@ -11,7 +11,7 @@ CXX = g++ -std=c++11
 INCLUDES = -I. -I/opt/local/include 
 
 # OPTFLAGS = -g
-OPTFLAGS = -O2 -g -pg
+OPTFLAGS = -O2 -g
 # OPTFLAGS = -O2 -g -DNDEBUG
 
 LDFLAGS = -L/opt/local/lib
@@ -120,6 +120,6 @@ realclean: clean
 $(LIB_OBJS): $(LIB_HEADERS)
 $(ALGEBRA_OBJS): $(ALGEBRA_HEADERS) $(LIB_HEADERS)
 $(KNOTKIT_OBJS) main.o mpimain.o kk.o: $(KNOTKIT_HEADERS) $(ALGEBRA_HEADERS) $(LIB_HEADERS) $(PERIODICITY_HEADERS)
-$(PERIODICITY_OBJS) : $(PERIODICITY_HEADERS)
+$(PERIODICITY_OBJS) : $(PERIODICITY_HEADERS) $(ALGEBRA_HEADERS) $(LIB_HEADERS)
 
 mpimain.o mpi_aux.o: mpi_aux.h

algebra/Z.h

@@ -14,6 +14,7 @@ class Z
 
  private:  
   std::shared_ptr<mpz_class> impl;
+#ifdef DEBUG_Z
   void write_state() const {
     std::cout << "I store the following value " << *this << "\n";
     std::cout << "Number of objects pointing to the same value " << impl.use_count() << "\n";
@@ -22,6 +23,7 @@ class Z
     /* std::cout << "Size of std::shared_ptr<mpz_class> " << sizeof(impl) << "\n"; */
     /* std::cout << "Size of mpz_class " << sizeof(*impl) << "\n"; */
   }
+#endif
   
  public:
   Z() : impl(new mpz_class) {
@@ -52,6 +54,12 @@ class Z
 #ifdef DEBUG_Z
     std::cout << "Z(const Z& z)" << "\n";
     write_state();
+#endif
+  }
+  Z(copy, const Z& z) : impl(z.impl) {
+#ifdef DEBUG_Z
+    std::cout << "Z(COPY, const Z& z)\n";
+    write_state();
 #endif
   }
   Z(Z&& z) : impl(std::move(z.impl)) {
@@ -115,6 +123,17 @@ class Z
     return *impl.get() < *z.impl.get();
   }
 
+  bool operator <= (const Z& z) const {
+    return *this < z || *this == z;
+  }
+
+  bool operator > (const Z& z) const {
+    return !(*this <= z);
+  }
+  bool operator >= (const Z& z) const {
+    return !(*this < z);
+  }
+
   bool is_unit () const
   {
     return *this == 1 || *this == -1;

algebra/multivariate_laurentpoly.h

@@ -548,7 +548,6 @@ multivariate_laurentpoly<T> invert_variable(const multivariate_laurentpoly<T>& p
   return result;
 }
 
-<<<<<<< HEAD
 template<class T>
 multivariate_laurentpoly<T>
 multivariate_laurentpoly<T>::evaluate(T val, unsigned index) const {
@@ -569,6 +568,4 @@ multivariate_laurentpoly<T>::evaluate(T val, unsigned index) const {
   return res;
 }
 
-=======
->>>>>>> 297e02dc05d79c4546d04aeb720680a9c231284e
 #endif // _KNOTKIT_ALGEBRA_MULTIVARIATE_LAURENPOLY_H

kk.cpp

@@ -1,10 +1,6 @@
-
 #include <knotkit.h>
-<<<<<<< HEAD
 #include <periodicity.h>
 #include <fstream>
-=======
->>>>>>> 297e02dc05d79c4546d04aeb720680a9c231284e
 #include <vector>
 #include <utility>
 
@@ -322,6 +318,17 @@ int compute_s_inv(knot_diagram& kd) {
 }
 
 void check_periodicity(std::string out_file) {
+  if(periodicity_test == "all") {
+    Kh_periodicity_checker Kh_pc(kd);
+    Przytycki_periodicity_checker P_pc(Kh_pc.get_KhP().evaluate(-1, eval_index));
+    for(auto& p : primes_list) {
+      std::cout << "Przytycki's criterion: "
+		<< P_pc(p) << std::endl
+		<< "Kh criterion: "
+		<< Kh_pc(p) << std::endl;
+    }
+  }
+  else {
   if(period == 2 || period == 3) {
     std::cout << "Sorry, the criteria don't work for period "
 	      << period << "...\n";
@@ -332,75 +339,28 @@ void check_periodicity(std::string out_file) {
     std::cout << "For now you can only check periodicity for primes up to 31..." << "\n";
     exit(EXIT_FAILURE);
   }
-  std::ofstream out(out_file);
+  std::ofstream out(out_file, std::ios_base::app);
   
-  if(periodicity_test == "all") {
+  if(periodicity_test == "Przytycki") {
-    std::cout << "I will perform both test..." << "\n";
+    Przytycki_periodicity_checker P_pc(compute_jones(kd));
-  }
-  else if(periodicity_test == "Przytycki") {
-    switch(period) {
-    case 5: {
-      Przytycki_periodicity_checker<5> pcc;
     if(out_file.size() != 0)
-	out << pcc(compute_jones(kd)) << std::endl;
+      out << P_pc(period) << std::endl;
     else
-	std::cout << pcc(compute_jones(kd)) << std::endl;
+      std::cout << P_pc(period) << std::endl;
-      break;
-    }
-    case 7: {
-      Przytycki_periodicity_checker<7> pcc;
-      if(out_file.size() != 0)
-	out << pcc(compute_jones(kd)) << std::endl;
-      else
-	std::cout << pcc(compute_jones(kd)) << std::endl;
-      break;
-    }
-    case 11: {
-      Przytycki_periodicity_checker<11> pcc;
-      if(out_file.size() != 0)
-	out << pcc(compute_jones(kd)) << std::endl;
-      else
-	std::cout << pcc(compute_jones(kd)) << std::endl;
-      break;
-    }
-    case 13: {
-      Przytycki_periodicity_checker<13> pcc;
-      if(out_file.size() != 0)
-	out << pcc(compute_jones(kd)) << std::endl;
-      else
-	std::cout << pcc(compute_jones(kd)) << std::endl;
-      break;
-    }
-    case 17: {
-      Przytycki_periodicity_checker<17> pcc;
-      if(out_file.size() != 0)
-	out << pcc(compute_jones(kd)) << std::endl;
-      else
-	std::cout << pcc(compute_jones(kd)) << std::endl;
-      break;
-    }
-    case 19: {
-      Przytycki_periodicity_checker<19> pcc;
-      if(out_file.size() != 0)
-	out << pcc(compute_jones(kd)) << std::endl;
-      else
-	std::cout << pcc(compute_jones(kd)) << std::endl;
-      break;
-    }
-    }
   }
   else if(periodicity_test == "Kh") {
-    Kh_periodicity_checker pc(kd);
+    Kh_periodicity_checker Kh_pc(kd);
     if(out_file.size() != 0)
-      out << pc(period) << std::endl;
+      out << Kh_pc(period) << std::endl;
     else
-      std::cout << pc(period) << std::endl;
+      std::cout << Kh_pc(period) << std::endl;
   }
   else {
       std::cout << "Sorry, I don't recognize this option..." << "\n";
       exit(EXIT_FAILURE);
     }
   }
+}
 
 
 template<class R> void
@@ -726,7 +686,7 @@ main (int argc, char **argv)
       else if(!strcmp (argv[i], "-p")) {
 	i++;
 	if(i == argc) {
-	  fprintf (stderr, "error: missing argument to option `-o'\n");
+	  fprintf (stderr, "error: missing argument to option `-p'\n");
 	  exit (EXIT_FAILURE);
 	}
 	period = std::stoi(argv[i]);
@@ -734,7 +694,7 @@ main (int argc, char **argv)
       else if (!strcmp(argv[i], "-t")) {
 	i++;
 	if(i == argc) {
-	  fprintf (stderr, "error: missing argument to option `-o'\n");
+	  fprintf (stderr, "error: missing argument to option `-t'\n");
 	  exit (EXIT_FAILURE);
 	}
 	periodicity_test = argv[i];
@@ -844,25 +804,6 @@ main (int argc, char **argv)
 	      << ") of " << knot <<  " = " << std::endl
 	      << khp << std::endl;
   }
-  else if(!strcmp(invariant, "periodicity_congruence")) {
-    if(!strcmp(field, "Z2")) {
-      // first we check whether period is a prime
-      if(period == 2 || period == 3) {
-	std::cerr << "The criterion does not work for period = " << period << "\n";
-	exit(EXIT_FAILURE);
-      }
-      auto result = find(primes.begin(), primes.end(), period);
-      if(result == primes.end()) {
-	std::cerr << "For now it is possible to check periodicity for primes up to 31" << "\n";
-	exit(EXIT_FAILURE);
-      }
-      check_periodicity_criterion<Z2>(kd,period);
-    }
-    else {
-      std::cerr << "error: for now this function is only defined for Z2 coefficients..." << "\n";
-      exit(EXIT_FAILURE);
-    }
-  }
   else
     {
       if (!strcmp (field, "Z2"))

lib/map.h

@@ -1,4 +1,3 @@
-
 /* wrapper for std::map */
 
 template<class K, class V>

periodicity.cpp

@@ -1,6 +1,150 @@
 #include <periodicity.h>
 #include <simplify_chain_complex.h>
 
+bool Przytycki_periodicity_checker::check(int period) const {
+  switch(period) {
+  case 5: {
+    periodic_congruence_checker<Zp<5>> pcc(5);
+    return pcc(jones_pol);
+  }
+  case 7: {
+    periodic_congruence_checker<Zp<7>> pcc(7);
+    return pcc(jones_pol);
+  }
+  case 11: {
+    periodic_congruence_checker<Zp<11>> pcc(11);
+    return pcc(jones_pol);
+  }
+  case 13: {
+    periodic_congruence_checker<Zp<13>> pcc(13);
+    return pcc(jones_pol);
+  }
+  case 17: {
+    periodic_congruence_checker<Zp<17>> pcc(17);
+    return pcc(jones_pol);
+  }
+  case 19: {
+    periodic_congruence_checker<Zp<19>> pcc(19);
+    return pcc(jones_pol);
+  }
+  }
+}
+
+std::string Przytycki_periodicity_checker::operator () (int period) const {
+  std::ostringstream res;
+  res << knot << ": period = " << period << ": "
+      << (check(period) ? "Maybe" : "No");
+  return res.str();
+}
+
+template<class T>
+polynomial_iterator<T>::polynomial_iterator(const multivariate_laurentpoly<T>& pol,
+					 start_pos sp) {
+  for(typename map<monomial, T>::const_iter i = pol.coeffs; i; ++i) {
+    monomials.push_back(i.key());
+    bounds.push_back(i.val());
+    current_pos.push_back(Z(0));
+  }
+
+  if(sp == start_pos::begin) {
+    level = 0;
+  }
+  else {
+    level = bounds.size();
+  }
+#ifndef NDEBUG
+  check_current_pos();
+#endif
+}
+
+#ifndef NDEBUG
+template<class T>
+void polynomial_iterator<T>::check_current_pos() {
+  assert(bounds.size() == monomials.size());
+  assert(bounds.size() == current_pos.size());
+  assert(level <= current_pos.size());
+  for(unsigned i = 0; i < current_pos.size(); i++) {
+    if(i < level) {
+      assert((current_pos[i] <= bounds[i]) &&
+	     Z(0) <= (current_pos[i]));
+    }
+    else if(i == level) {
+      if(level > 0)
+	assert(current_pos[i] <= bounds[i] && Z(0) < current_pos[i]);
+      else
+	assert((current_pos[i] <= bounds[i]) && Z(0) <= (current_pos[i]));	
+    }
+    else
+      assert(current_pos[i] == Z(0));
+  }
+}
+#endif // NDEBUG
+
+template<class T>
+polynomial_iterator<T>&
+polynomial_iterator<T>::operator ++ () {
+#ifndef NDEBUG
+  check_current_pos();
+#endif
+  if(level == monomials.size())
+    return *this;
+
+  unsigned i = 0;
+  while(i <= level) {
+#ifndef NDEBUG
+    check_current_pos();
+#endif
+    if(current_pos[i] < bounds[i]) {
+      current_pos[i] += 1;
+      break;
+    }
+    else {
+      if(i == level) {
+	if(level < monomials.size() - 1) {
+	  current_pos[i] = 0;
+	  current_pos[i+1] += 1;
+	  level++;
+	  break;
+	}
+	else {
+	  level++;
+	  break;
+	}
+      }
+      else {
+	current_pos[i] = 0;
+      }
+    }
+    i++;
+  }
+  return *this;
+}
+
+template<class T>
+multivariate_laurentpoly<T> polynomial_iterator<T>::operator *() const {
+  polynomial res;
+  for(unsigned i = 0; i <= level; i++) {
+    res += polynomial(current_pos[i], monomials[i]);
+  }
+  return res;
+}
+
+template<class T>
+std::string polynomial_iterator<T>::write_self() const {
+  std::ostringstream res;
+  res << "level = " << level << std::endl
+      << "monomials:" << std::endl;
+  for(auto& mon : monomials)
+    res << mon << std::endl;
+  res << "bounds: " << std::endl;
+  for(auto& b : bounds)
+    res << b << std::endl;
+  res << "current_pos: " << std::endl;
+  for(auto& pos : current_pos)
+    res << pos << std::endl;
+  return res.str();
+}
+
 void Kh_periodicity_checker::compute_knot_polynomials(knot_diagram& kd) {
 
   unsigned m = kd.num_components ();
@@ -57,8 +201,14 @@ void Kh_periodicity_checker::compute_quot() {
 	quot += polynomial(m1.second, m1.first);
 	polynomial p = polynomial(m1.second, m1.first) * mul;
 	diff -= p;
+	if(diff != 0)
 	  m1 = diff.tail();
+	else break;
       }
+      if(diff != 0)
+	m = diff.head();
+      else
+	break;
     }
     quot += polynomial(m.second, m.first);
     polynomial p = polynomial(m.second, m.first) * mul;
@@ -91,33 +241,6 @@ Kh_periodicity_checker::compute_quotient_and_remainder(const polynomial& quot,
   return std::make_pair(quotient, remainder);
 }
 
-std::list<multivariate_laurentpoly<Z>>
-Kh_periodicity_checker::generate_candidates(const polynomial &q) const {
-  std::list<polynomial> result;
-  Z size = 0;
-  map<monomial, Z>::const_iter i = q.coeffs;
-
-  for(int j = 0; Z(j) <= i.val(); j++) {
-    result.push_back(polynomial(Z(j), i.key()));
-    size += 1;
-  }
-
-  i++;
-  
-  for( ; i; i++) {
-    for(int j = 1; Z(j) <= i.val(); j++) {
-      std::list<polynomial> temp_list;
-      for_each(result.begin(), result.end(),
-	       [&result, &temp_list, j, i](const polynomial& p){
-		 temp_list.push_back(p + polynomial(Z(j), i.key()));
-	       });
-      result.splice(result.end(), temp_list);
-      size += 1;
-    }
-  }
-  return result;
-}
-
 bool Kh_periodicity_checker::check(const polynomial& q,
 				   const polynomial& r,
 				   int period) const {
@@ -126,26 +249,34 @@ bool Kh_periodicity_checker::check(const polynomial& q,
   if(q == 0) {
     return pcc(t.evaluate(-1,1));
   }
-  else {
-    // generate all polynomials
-    if(verbose) {
-      std::cout << "Generating all candidates..." << std::endl;
-    }
-    std::list<polynomial> candidates = generate_candidates(q);
-    // Z i = 0;
-    // for(auto& l : candidates) {
-    //   i += 1;
-    //   std::cout << i << ": " << l << std::endl;
-    // }
-    // and check each one of them
   if(verbose)
-      std::cout << "Checking congruences..." << std::endl;
+    std::cout << "Checking congruences..."; 
-    polynomial m = mul;
+  polynomial_iterator<Z> pi(q);
-    return any_of(candidates.begin(), candidates.end(),
+  polynomial_iterator<Z> pi_end(q, polynomial_iterator<Z>::start_pos::end);
-		  [&pcc, &t, &m, &period](const polynomial& p){
+  Z count = pi.get_count();
-		    return pcc((t + Z(period - 1) * m * p).evaluate(-1,1));
+  if(verbose)
-      });
+    std::cout << count << " candidates..." << std::endl;
+  Z step = 1;
+  if(count >= 32)
+    step = count / 32;
+  Z c = 0;
+  while(pi != pi_end) {
+    //std::cout << "pi: " << std::endl << pi;
+    polynomial temp = t + polynomial(period - 1) * mul * (*pi);
+    if(pcc(temp.evaluate(-1,1))) {
+      std::cout << "Candidates:" << std::endl
+  		<< "EKhP_1 = " << temp << std::endl
+  		<< "EKhP_" << (period - 1) << " = "
+  		<< (polynomial(period - 1) * mul *(r - *pi))
+  		<< std::endl;
+      return true;
     }
+    ++pi;
+    c += 1;
+    if(verbose && c % step == 0)
+      std::cout << c << "/" << count << "..." << std::endl;
+  }
+  return false;
 }
 
 std::string Kh_periodicity_checker::operator () (int period) const {

periodicity.h

@@ -5,9 +5,6 @@
 #include <sstream>
 #include <string>
 #include <vector>
-#include <utility>
-#include <list>
-#include <algorithm>
 
 extern bool verbose;
 extern const char* knot;
@@ -16,15 +13,12 @@ extern std::string periodicity_test;
 
 const std::vector<int> primes_list = {5, 7, 11, 13, 17, 19};
 
-enum class Test_type {
+const unsigned eval_index = 1;
-  Przytycki_criterion,
+const unsigned invert_index = 2;
-  BKP_criterion,
-  all
-};
 
 template<class T>
 class periodic_congruence_checker {
-public:
+protected:
   using polynomial = multivariate_laurentpoly<T>;
   using monomial = multivariate_laurent_monomial;
   
@@ -32,19 +26,20 @@ public:
   unsigned index;
 
   polynomial prepare_polynomial(const polynomial& pol) const {
-    return (pol - invert_variable(pol, index));
+    polynomial inv = invert_variable(pol, index);
+    return pol - inv;
   }
 
   bool reduce(const polynomial& pol) const;
 
 public:
   periodic_congruence_checker(int pprime = 5,
-			      unsigned ind = 2) :
+			      unsigned ind = invert_index) :
     prime(pprime),
     index(ind)
     {}
 
-  ~periodic_congruence_checker() {};
+  virtual ~periodic_congruence_checker() {};
 
   bool operator() (const polynomial& pol) const {
     return reduce(prepare_polynomial(pol));
@@ -61,32 +56,90 @@ bool periodic_congruence_checker<T>::reduce(const multivariate_laurentpoly<T>& p
     monomial mon = monomial(VARIABLE, index, c);
     res += polynomial(i.val(), mon);
   }
-  if(verbose)
+  // if(verbose)
-    std::cout << "reduced = " << res << "\n";
+  //   std::cout << "reduced = " << res << "\n";
   return res == 0;
 }
 
-template<unsigned p>
 class Przytycki_periodicity_checker {
-  using polynomial = multivariate_laurentpoly<Zp<p>>;
+  using polynomial = multivariate_laurentpoly<Z>;
   using monomial = multivariate_laurent_monomial;
 
-  periodic_congruence_checker<Zp<p>> cong_checker;
+  polynomial jones_pol;
+
+  bool check(int period) const;
 
  public:
-  Przytycki_periodicity_checker() : cong_checker(p) {}
+  Przytycki_periodicity_checker(polynomial j) : jones_pol(j) {}
     
   ~Przytycki_periodicity_checker() {}
   
-  std::string operator() (const polynomial& pol) const {
+  std::string operator() (int period) const;
-    std::ostringstream out;
+};
-    out << knot << ": period = " << p
+
-	<< ": "
+template<class T>
-	<< (cong_checker(pol) ? "Maybe" : "No");
+class polynomial_iterator {
-    return out.str();
+  using polynomial = multivariate_laurentpoly<T>;
+  using monomial = multivariate_laurent_monomial;
+
+  std::vector<monomial> monomials;
+  std::vector<T> bounds;
+  std::vector<T> current_pos;
+  unsigned level;
+
+  void check_current_pos();
+
+public:
+  enum class start_pos { begin, end };
+  
+  polynomial_iterator(const polynomial& init,
+		      start_pos sp = start_pos::begin);
+  polynomial_iterator(const polynomial_iterator& pi) =default;
+  polynomial_iterator(polynomial_iterator&& pi) =default;
+
+  ~polynomial_iterator() {}
+
+  polynomial_iterator& operator = (const polynomial_iterator& pi) =default;
+  polynomial_iterator& operator = (polynomial_iterator&& pi) =default;
+
+  polynomial_iterator& operator ++();
+
+  bool operator == (const polynomial_iterator& pi) const {
+    if(level == monomials.size() || pi.level == pi.monomials.size()) {
+      return level == pi.level &&
+	monomials == pi.monomials &&
+	bounds == pi.bounds;
+    }
+    else {
+      return level == pi.level &&
+	bounds == pi.bounds &&
+	monomials == pi.monomials &&
+	current_pos == pi.current_pos;
+    }
+  }
+  bool operator != (const polynomial_iterator& pi) const {
+    return !(*this == pi);
+  }
+  polynomial operator*() const;
+
+  T get_count() const {
+    Z res = 1;
+    for(auto& v : bounds)
+      res *= (v + 1);
+    return res;
+  }
+
+  std::string write_self() const;
+  friend inline std::ostream& operator << (std::ostream& os, const polynomial_iterator& pi) {
+    return os << pi.write_self();
   }
 };
 
+template<class T>
+std::ostream& operator << (std::ostream& os, const polynomial_iterator<T>& pi) {
+  return os << *pi;
+}
+
 class Kh_periodicity_checker {
   using polynomial = multivariate_laurentpoly<Z>;
   using monomial = multivariate_laurent_monomial;
@@ -101,7 +154,7 @@ class Kh_periodicity_checker {
   void compute_quot();
   std::pair<polynomial, polynomial> compute_quotient_and_remainder(const polynomial& p,
 								   int period) const;
-  std::list<polynomial> generate_candidates(const polynomial& q) const;
+  // std::list<polynomial> generate_candidates(const polynomial& q) const;
   bool check(const polynomial& q, const polynomial& r, int period) const;
 
  public:
@@ -119,6 +172,14 @@ class Kh_periodicity_checker {
   ~Kh_periodicity_checker() {}
 
   std::string operator () (int period) const;
+
+  polynomial get_KhP() const {
+    return khp;
+  }
+
+  polynomial get_LeeP() const {
+    return leep;
+  }
 };
 
 #endif // _KNOTKIT_PERIODICITY_H