diff --git a/.gitignore b/.gitignore
index 093d3f1..62f5ee3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -11,3 +11,4 @@ GTAGS
 /parallel.cmd.[eo]*
 /save
 */save
+/org
diff --git a/Makefile b/Makefile
index 18109a1..3bb09c6 100644
--- a/Makefile
+++ b/Makefile
@@ -27,7 +27,9 @@ KNOTKIT_OBJS = planar_diagram.o dt_code.o knot_diagram.o cube.o steenrod_square.
   knot_parser/knot_parser.o knot_parser/knot_scanner.o \
   rd_parser/rd_parser.o rd_parser/rd_scanner.o
 
-COMMON_OBJS = $(KNOTKIT_OBJS) $(ALGEBRA_OBJS) $(LIB_OBJS) 
+PERIODICITY_OBJS = periodicity.o
+
+COMMON_OBJS = $(KNOTKIT_OBJS) $(ALGEBRA_OBJS) $(LIB_OBJS) $(PERIODICITY_OBJS)
 
 LIB_HEADERS = lib/lib.h lib/show.h lib/refcount.h lib/pair.h lib/maybe.h lib/vector.h \
   lib/set_wrapper.h lib/set.h lib/hashset.h \
@@ -44,6 +46,8 @@ KNOTKIT_HEADERS = knotkit.h planar_diagram.h dt_code.h knot_diagram.h \
   smoothing.h cobordism.h cube.h steenrod_square.h \
   spanning_tree_complex.h cube_impl.h sseq.h simplify_chain_complex.h
 
+PERIODICITY_HEADERS = periodicity.h
+
 LIBS = -lgmpxx -lgmp  -lz
 
 all: kk
@@ -55,7 +59,7 @@ all: kk
 	$(CXX) -c $(CXXFLAGS) $< -o $@
 
 kk: kk.o $(COMMON_OBJS)
-	$(CXX) $(LDFLAGS) -o kk $^ $(LIBS)
+	$(CXX) $(LDFLAGS) -pg -o kk $^ $(LIBS)
 
 main: main.o $(COMMON_OBJS)
 	$(CXX) $(LDFLAGS) -o main $^ $(LIBS)
@@ -115,6 +119,7 @@ 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)
+$(KNOTKIT_OBJS) main.o mpimain.o kk.o: $(KNOTKIT_HEADERS) $(ALGEBRA_HEADERS) $(LIB_HEADERS) $(PERIODICITY_HEADERS)
+$(PERIODICITY_OBJS) : $(PERIODICITY_HEADERS)
 
 mpimain.o mpi_aux.o: mpi_aux.h
diff --git a/algebra/Q.h b/algebra/Q.h
index 24c8509..b8f2510 100644
--- a/algebra/Q.h
+++ b/algebra/Q.h
@@ -13,47 +13,8 @@ class Q
   using linear_combination_const_iter = ::linear_combination_const_iter<Q>;
 
  private:
-  // enum steal { STEAL };
-  // enum copy { COPY };
-  // class Q_impl : public refcounted
-  // {
-  // public:
-  //   mpq_t x;
-    
-  // public:
-  //   Q_impl () { mpq_init (x); }
-  //   Q_impl (int init)
-  //   {
-  //     mpq_init (x);
-  //     mpq_set_si (x, init, 1);
-  //   }
-    
-  //   Q_impl (copy, mpq_srcptr init)
-  //   {
-  //     mpq_init (x);
-  //     mpq_set (x, init);
-  //   }
-    
-  //   Q_impl (steal, mpq_srcptr init) { x[0] = *init; }
-  //   Q_impl (reader &r)
-  //   {
-  //     mpq_init (x);
-  //     r.read_mpz (mpq_numref (x));
-  //     r.read_mpz (mpq_denref (x));
-  //   }
-    
-  //   ~Q_impl () { mpq_clear (x); }
-    
-  //   void write_self (writer &w) const
-  //   {
-  //     w.write_mpz (mpq_numref (x));
-  //     w.write_mpz (mpq_denref (x));
-  //   }
-  // };
-  
-  // ptr<Q_impl> impl;
   std::shared_ptr<mpq_class> impl;
-    void write_state() const {
+  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";
     /* std::cout << "I point to " << impl.get() << "\n"; */
@@ -246,6 +207,11 @@ class Q
   int get_count() const {
     return impl.use_count();
   }
+  Z get_num() const {
+    return Z(impl.get()->get_num());
+  }
+  Z get_den() const {
+    return Z(impl.get()->get_den());
+  }
 };
-
 #endif // _KNOTKIT_ALGEBRA_Q_H
diff --git a/algebra/Z.h b/algebra/Z.h
index b15856a..db61310 100644
--- a/algebra/Z.h
+++ b/algebra/Z.h
@@ -145,6 +145,24 @@ class Z
     }
   }
 
+  Z operator % (const Z& z) const {
+    if(Z(0) < z) {
+      // return Z(*impl.get() % *z.impl.get());
+      mpz_class r;
+      mpz_fdiv_r(r.get_mpz_t(), impl.get()->get_mpz_t(), z.impl.get()->get_mpz_t());
+      return Z(r);
+    }    
+    else
+      return *this;
+  }
+
+  Z operator % (unsigned p) const {
+    if(p != 0)
+      return operator % (Z(p));
+    else
+      return *this;
+  }
+
   Z recip () const {
     assert(is_unit());
     return *this;
@@ -234,6 +252,9 @@ class Z
   int get_count() const {
     return impl.use_count();
   }
-};
 
+  unsigned get_ui() const {
+    return impl.get()->get_ui();
+  }
+};
 #endif // _KNOTKIT_ALGEBRA_Z_H
diff --git a/algebra/Z2.h b/algebra/Z2.h
index 8cb5df5..c5c2ba6 100644
--- a/algebra/Z2.h
+++ b/algebra/Z2.h
@@ -1,5 +1,6 @@
 #ifndef KNOTKIT_ALGEBRA_Z2_H
 #define KNOTKIT_ALGEBRA_Z2_H
+#include <algebra/Z.h>
 #include <iostream>
 #include <tuple>
 #include <cassert>
@@ -8,7 +9,6 @@ class Z2
 {
  public:
   using linear_combination = ::linear_combination<Z2>;
-  // typedef linear_combination_iter<Z2> linear_combination_iter;
   using linear_combination_const_iter = ::linear_combination_const_iter<Z2>;
   
  private:
@@ -24,6 +24,7 @@ class Z2
     x.v = 0;
   }
   Z2 (reader &r) { v = r.read_bool (); }
+  Z2 (const Z& z) : v((z % 2).get_ui()) {}
   ~Z2 () { }
   
   Z2& operator = (const Z2& x) { v = x.v; return *this; }
@@ -88,5 +89,4 @@ class Z2
   void show_self () const { std::cout << *this; }
   void display_self () const { std::cout << *this << "\n"; }
 };
-
 #endif //KNOTKIT_ALGEBRA_Z2_H
diff --git a/algebra/Zp.h b/algebra/Zp.h
index 80414a0..3de8b6d 100644
--- a/algebra/Zp.h
+++ b/algebra/Zp.h
@@ -1,6 +1,6 @@
 #ifndef KNOTKIT_ALGEBRA_ZP_H
 #define KNOTKIT_ALGEBRA_ZP_H
-
+#include <algebra/Z.h>
 #include <iostream>
 #include <tuple>
 #include <cassert>
@@ -32,6 +32,7 @@ class Zp
   Zp (Zp&& x) : v(std::move(x.v)) {
     x.v = 0;
   }
+  Zp (const Z& z) : v((z % p).get_ui()) {}
   ~Zp () { }
   
   Zp& operator = (const Zp& x) { v = x.v; return *this; }
@@ -123,5 +124,4 @@ class Zp
   void show_self () const { std::cout << v << "(" << p << ")"; }
   void display_self () const { std::cout << *this << "\n"; }
 };
-
 #endif //KNOTKIT_ALGEBRA_ZP_H
diff --git a/algebra/algebra.h b/algebra/algebra.h
index b23bd79..325e0ba 100644
--- a/algebra/algebra.h
+++ b/algebra/algebra.h
@@ -1,3 +1,5 @@
+#ifndef _KNOTKIT_ALGEBRA_H
+#define _KNOTKIT_ALGEBRA_H
 #include <lib/lib.h>
 
 #ifdef DEBUG_ALGEBRA
@@ -73,3 +75,5 @@ enum variable { VARIABLE };
 
 #include <algebra/module.h>
 #include <algebra/linear_combination.h>
+
+#endif // _KNOTKIT_ALGEBRA_H
diff --git a/algebra/fraction_field.h b/algebra/fraction_field.h
index 3f648b6..ea3a6d2 100644
--- a/algebra/fraction_field.h
+++ b/algebra/fraction_field.h
@@ -1,3 +1,5 @@
+#ifndef _KNOTKIT_ALGEBRA_FRACTION_FIELD_H
+#define _KNOTKIT_ALGEBRA_FRACTION_FIELD_H
 
 template<class T> class fraction_field
 {
@@ -243,3 +245,5 @@ fraction_field<T>::show_self () const
       printf (")");
     }
 }
+
+#endif // _KNOTKIT_ALGEBRA_FRACTION_FIELD_H
diff --git a/algebra/grading.h b/algebra/grading.h
index cbd3972..4310657 100644
--- a/algebra/grading.h
+++ b/algebra/grading.h
@@ -1,4 +1,5 @@
-
+#ifndef _KNOTKIT_ALGEBRA_GRADINGS_H
+#define _KNOTKIT_ALGEBRA_GRADINGS_H
 class grading
 {
  public:
@@ -8,6 +9,9 @@ class grading
   grading () : h(0), q(0) { }
   grading (int h_, int q_) : h(h_), q(q_) { }
   grading (const grading &gr) : h(gr.h), q(gr.q) { }
+  grading (grading&& gr) : h(std::move(gr.h)), q(std::move(gr.q)) {
+    gr.h = gr.q = 0;
+  }
   grading (reader &r)
   {
     h = r.read_int ();
@@ -17,6 +21,12 @@ class grading
   ~grading () { }
   
   grading &operator = (const grading &gr) { h = gr.h; q = gr.q; return *this; }
+  grading& operator = (grading&& gr) {
+    h = std::move(gr.h);
+    q = std::move(gr.q);
+    gr.h = gr.q = 0;
+    return *this;
+  }
   
   grading operator + (const grading &gr) const
   {
@@ -61,3 +71,5 @@ class grading
   void show_self () const;
   void display_self () const;
 };
+
+#endif // _KNOTKIT_ALGEBRA_GRADINGS_H
diff --git a/algebra/linear_combination.h b/algebra/linear_combination.h
index 373a0b8..b9af52b 100644
--- a/algebra/linear_combination.h
+++ b/algebra/linear_combination.h
@@ -1,4 +1,5 @@
-
+#ifndef _KNOTKIT_ALGEBRA_LINEAR_COMBINATIONS_H
+#define _KNOTKIT_ALGEBRA_LINEAR_COMBINATIONS_H
 template<class R> class linear_combination_const_iter;
 
 template<class R>
@@ -512,3 +513,5 @@ class linear_combination_const_iter<Z2>
   unsigned key () const { return i.val (); }
   Z2 val () { return Z2 (1); }
 };
+
+#endif // _KNOTKIT_ALGEBRA_LINEAR_COMBINATIONS_H
diff --git a/algebra/module.h b/algebra/module.h
index acd217b..e824ab9 100644
--- a/algebra/module.h
+++ b/algebra/module.h
@@ -1,4 +1,5 @@
-
+#ifndef _KNOTKIT_ALGEBRA_MODULE_H
+#define _KNOTKIT_ALGEBRA_MODULE_H
 template<class R> class mod_map;
 
 template<class R> class mod_span;
@@ -85,8 +86,10 @@ class module : public refcounted
   ptr<const quotient_module<R> > quotient (ptr<const free_submodule<R> > m) const;
   
   ptr<const free_submodule<R> > submodule (const mod_span<R> &span) const;
-  
+
   multivariate_laurentpoly<Z> free_poincare_polynomial () const;
+  template<unsigned p>
+  multivariate_laurentpoly<Zp<p>> free_poincare_polynomial () const;
   multivariate_laurentpoly<Z> free_delta_poincare_polynomial () const;
   multivariate_laurentpoly<Z> free_ell_poincare_polynomial () const;
   
@@ -1476,6 +1479,20 @@ module<R>::free_poincare_polynomial () const
   return r;
 }
 
+template<class R> template <unsigned p> multivariate_laurentpoly<Zp<p>>
+module<R>::free_poincare_polynomial() const {
+  multivariate_laurentpoly<Zp<p>> r;
+  for (unsigned i = 1; i <= free_rank (); i ++)
+    {
+      grading hq = generator_grading (i);
+      multivariate_laurent_monomial m;
+      m.push_exponent (1, hq.h);
+      m.push_exponent (2, hq.q);
+      r.muladdeq (1, m);
+    }
+  return r;
+}
+
 template<class R> multivariate_laurentpoly<Z>
 module<R>::free_delta_poincare_polynomial () const
 {
@@ -2130,3 +2147,5 @@ reader::read_mod ()
       return module<R>::reader_id_module(id);
     }
 }
+
+#endif // _KNOTKIT_ALGEBRA_MODULE_H
diff --git a/algebra/multivariate_laurentpoly.h b/algebra/multivariate_laurentpoly.h
index a26122c..6c9d548 100644
--- a/algebra/multivariate_laurentpoly.h
+++ b/algebra/multivariate_laurentpoly.h
@@ -1,4 +1,7 @@
-
+#ifndef _KNOTKIT_ALGEBRA_MULTIVARIATE_LAURENPOLY_H
+#define _KNOTKIT_ALGEBRA_MULTIVARIATE_LAURENPOLY_H
+#include <sstream>
+#include <string>
 /* multivariate polynomial in a (vector) variable x with coefficients
    in T. */
 
@@ -133,17 +136,31 @@ class multivariate_laurent_monomial
     if (e != 0)
       m.push (j, e);
   }
-  
+
+  std::string to_string() const {
+    std::ostringstream res;
+    for (map<unsigned, int>::const_iter i = m; i; i ++) {
+      assert (i.val () != 0);
+      if (i.val () == 1) {
+	res << "x" << std::to_string(i.key()); 
+      }
+      else {
+	res << "x"
+	    << std::to_string(i.key())
+	    << "^"
+	    << std::to_string(i.val());
+      }
+    }
+    return res.str();
+  }
+
+  friend std::ostream& operator << (std::ostream& os, const multivariate_laurent_monomial& m) {
+    return os << m.to_string();
+  }
+    
   void show_self () const
   {
-    for (map<unsigned, int>::const_iter i = m; i; i ++)
-      {
-	assert (i.val () != 0);
-	if (i.val () == 1)
-	  printf ("x%d", i.key ());
-	else
-	  printf ("x%d^%d", i.key (), i.val ());
-      }
+    std::cout << *this;
   }
   
   void write_self (writer &w) const { write (w, m); }
@@ -161,12 +178,11 @@ template<class T>
 class multivariate_laurentpoly
 {
  public:
-  typedef ::linear_combination<multivariate_laurentpoly<T> > linear_combination;
-  typedef ::linear_combination_const_iter<multivariate_laurentpoly<T> >
-    linear_combination_const_iter;
+  using linear_combination = ::linear_combination<multivariate_laurentpoly<T>>;
+  using linear_combination_const_iter = ::linear_combination_const_iter<multivariate_laurentpoly<T> >;
   
  public:
-  typedef multivariate_laurent_monomial monomial;
+  using monomial = multivariate_laurent_monomial;
   
   map<monomial, T> coeffs;
   
@@ -216,13 +232,15 @@ class multivariate_laurentpoly
     if (c != 0)
       coeffs.push (m, c);
   }
+
+  template<class U>
+  multivariate_laurentpoly (const multivariate_laurentpoly<U>& pol);
   
   multivariate_laurentpoly (const multivariate_laurentpoly &p) : coeffs(p.coeffs) { }
   multivariate_laurentpoly (copy, const multivariate_laurentpoly &p)
     // ??? COPY2?
     : coeffs(COPY, p.coeffs)
-  {
-  }
+  {}
 
   multivariate_laurentpoly (reader &r) : coeffs(r) { }
   
@@ -294,6 +312,7 @@ class multivariate_laurentpoly
   
   unsigned card () const { return coeffs.card (); }
   pair<monomial, T> head () const { return coeffs.head (); }
+  pair<monomial, T> tail() const { return coeffs.tail(); }
   
   multivariate_laurentpoly &operator += (const multivariate_laurentpoly &p);
   multivariate_laurentpoly &operator -= (const multivariate_laurentpoly &p);
@@ -345,18 +364,33 @@ class multivariate_laurentpoly
 #ifndef NDEBUG
   void check () const;
 #endif
+
+  multivariate_laurentpoly evaluate(T val, unsigned index) const;
   
   static void show_ring ()
   {
     T::show_ring ();
     printf ("[x_1^+/-1, ..., x_n^+/-1]");
   }
-  
-  void display_self () const { show_self (); newline (); }
-  void show_self () const;
+
+  std::string to_string() const;
+  void display_self () const {
+    std::cout << *this << "\n";
+  }
+  void show_self () const {
+    std::cout << *this;
+  }
   void write_self (writer &w) const { write (w, coeffs); }
 };
 
+template<class T> template<class U>
+multivariate_laurentpoly<T>::multivariate_laurentpoly(const multivariate_laurentpoly<U>& pol) {
+  for(typename map<monomial, U>::const_iter i = pol.coeffs; i; i++) {
+    if(T(i.val()) != 0)
+      coeffs.push(i.key(), T(i.val()));
+  }
+}
+
 template<class T> multivariate_laurentpoly<T>
 operator * (const T &s, const multivariate_laurentpoly<T> &p)
 {
@@ -444,42 +478,50 @@ multivariate_laurentpoly<T>::check () const
 }
 #endif
 
-template<class T> void
-multivariate_laurentpoly<T>::show_self () const
-{
+template<class T>
+std::string multivariate_laurentpoly<T>::to_string() const {
+  std::ostringstream res;
   unsigned first = 1;
-  for (typename map<monomial, T>::const_iter i = coeffs; i; i ++)
-    {
-      monomial m = i.key ();
-      T c = i.val ();
+  for (typename map<monomial, T>::const_iter i = coeffs; i; i ++) {
+    monomial m = i.key ();
+    T c = i.val ();
+    assert (c != 0);
+
+    if (first)
+      first = 0;
+    else
+      res << " + ";
       
-      assert (c != 0);
-      
-      if (first)
-	first = 0;
+    if (m == 1) {
+      if (c == 1)
+	res << "1";
       else
-	printf (" + ");
-      
-      if (m == 1)
-	{
-	  if (c == 1)
-	    printf ("1");
-	  else
-	    show (c);
-	}
-      else
-	{
-	  if (c != 1)
-	    {
-	      show (c);
-	      printf ("*");
-	    }
-	  
-	  show (m);
-	}
+	res << c;
     }
+    else {
+      if (c != 1) {
+	res << c << "*";
+      }
+      res << m;
+    }
+  }
   if (first)
-    printf ("0");
+    res << "0";
+  return res.str();
+}
+
+template<class T>
+std::ostream& operator << (std::ostream& os, const multivariate_laurentpoly<T>& pol) {
+  return os << pol.to_string();
+}
+
+template<class T, class U>
+multivariate_laurentpoly<T> reduce(const multivariate_laurentpoly<U>& pol) {
+  multivariate_laurentpoly<T> res;
+  for(typename map<unsigned, U>::const_iter i = pol.coeffs; i; i++) {
+    res += multivariate_laurentpoly<T>(i.val(), i.key());
+  }
+  return res;
 }
 
 // functions below were added to verify several periodicity criteria
@@ -507,29 +549,23 @@ multivariate_laurentpoly<T> invert_variable(const multivariate_laurentpoly<T>& p
 }
 
 template<class T>
-bool check_przytycki_cong(const multivariate_laurentpoly<T>& pol, const int prime,
-			  const int exponent = 1, const unsigned index = 1) {
-  multivariate_laurentpoly<T> result;
-  for(typename map<multivariate_laurent_monomial, T>::const_iter i = pol.coeffs; i; i++) {
-    multivariate_laurent_monomial mon = i.key();
-    T c = i.val();
-    multivariate_laurent_monomial mon2;
-    for(typename map<unsigned, int>::const_iter i = mon.m; i; i++) {
-      // still need to handle the coefficient
-      // i.e. we need to take c % p
-      if(i.key() == index) {
-	int v = i.val() % (2 * prime);
-	if(v < 0) v += (2 * prime);
-	mon2 *= multivariate_laurent_monomial(VARIABLE, i.key(), v);
-	c.display_self();
-	printf("Old: key = %d, val = %d\n", i.key(), i.val());
-	printf("New: key = %d, val = %d\n", i.key(), v);
-      }
-      else
-	mon2 *= multivariate_laurent_monomial(VARIABLE, i.key(), i.val());
+multivariate_laurentpoly<T>
+multivariate_laurentpoly<T>::evaluate(T val, unsigned index) const {
+  using polynomial = multivariate_laurentpoly<T>;
+  using monomial = multivariate_laurent_monomial;
+  polynomial res;
+  for(typename map<monomial, T>::const_iter i = coeffs; i; i++) {
+    if(i.key().m % index) {
+      int exp = i.key().m[index];
+      monomial mon = i.key();
+      mon.m.yank(index);
+      polynomial temp = polynomial(i.val() * pow(val, exp), mon);
+      res += temp;
     }
-    result += multivariate_laurentpoly<T>(c, mon2);
+    else
+      res += polynomial(i.val(), i.key());
   }
-  result.display_self();
-  return true;
+  return res;
 }
+
+#endif // _KNOTKIT_ALGEBRA_MULTIVARIATE_LAURENPOLY_H
diff --git a/algebra/multivariate_polynomial.h b/algebra/multivariate_polynomial.h
index c35bd32..a30bc0d 100644
--- a/algebra/multivariate_polynomial.h
+++ b/algebra/multivariate_polynomial.h
@@ -1,4 +1,5 @@
-
+#ifndef _KNOTKIT_ALGEBRA_MULTIVARIATE_POLYNOMIAL_H
+#define _KNOTKIT_ALGEBRA_MULTIVARIATE_POLYNOMIAL_H
 /* multivariate polynomial in a (vector) variable x with coefficients
    in T. */
 
@@ -499,3 +500,5 @@ multivariate_polynomial<T, n>::show_self () const
   if (first)
     printf ("0");
 }
+
+#endif // _KNOTKIT_ALGEBRA_MULTIVARIATE_POLYNOMIAL_H
diff --git a/algebra/polynomial.h b/algebra/polynomial.h
index 4d1771f..c4bae40 100644
--- a/algebra/polynomial.h
+++ b/algebra/polynomial.h
@@ -1,4 +1,5 @@
-
+#ifndef _KNOTKIT_ALGEBRA_POLYNOMIAL_H
+#define _KNOTKIT_ALGEBRA_POLYNOMIAL_H
 /* univariate polynomial in a single variable `x' with coefficients in
    T. */
 
@@ -512,3 +513,5 @@ class polynomial<Z2>
   void display_self () const { show_self (); newline (); }
   void show_self () const;
 };
+
+#endif // _KNOTKIT_ALGEBRA_POLYNOMIAL_H
diff --git a/kk.cpp b/kk.cpp
index 2ddb2a7..c867d9c 100644
--- a/kk.cpp
+++ b/kk.cpp
@@ -1,4 +1,8 @@
 #include <knotkit.h>
+#include <periodicity.h>
+#include <fstream>
+#include <vector>
+#include <utility>
 
 const char *program_name;
 
@@ -70,9 +74,10 @@ void tex_footer ()
 const char *knot = 0;
 const char *invariant = 0;
 const char *field = "Z2";
+std::string periodicity_test = "Przytycki"; 
+int period = 5;
 
 knot_diagram kd;
-
 bool reduced = 0;
 
 class hg_grading_mapper
@@ -250,45 +255,197 @@ compute_gss ()
   tex_footer ();
 }
 
-multivariate_laurentpoly<Z> compute_jones(const knot_diagram& k, bool reduced) {
-  cube<Z2> c(kd, reduced);
-  ptr<const module<Z2> > C = c.khC;
-  mod_map<Z2> d = c.compute_d(1, 0, 0, 0, 0);
-  chain_complex_simplifier<Z2> s(C, d, maybe<int>(1), maybe<int>(0));
+template<class R>
+multivariate_laurentpoly<Z> compute_khp(const knot_diagram& k, bool reduced = false) {
+  cube<R> c (kd, reduced);
+  ptr<const module<R> > C = c.khC;
+  mod_map<R> d = c.compute_d (1, 0, 0, 0, 0);
+    
+  unsigned m = kd.num_components ();
+  hg_grading_mapper mapper (m);
+      
+  chain_complex_simplifier<R> s (C, d,
+				     maybe<int> (1), maybe<int> (0));
   C = s.new_C;
   d = s.new_d;
-  multivariate_laurentpoly<Z> jones;
-  for(uint i = 1; i <= c.n_generators; ++i) {
-    grading gr = c.compute_generator_grading(i);
-    if(gr.h % 2 == 0) {
-      jones += multivariate_laurentpoly<Z>(1, VARIABLE, 1, gr.q);
-    }
-    else {
-      jones += multivariate_laurentpoly<Z>(-1, VARIABLE, 1, gr.q);
-    }
-  }
-  return jones;
+  return C->free_poincare_polynomial();
+}
+
+multivariate_laurentpoly<Z> compute_jones(const knot_diagram& k, bool reduced = false) {
+  return compute_khp<Z2>(k, reduced).evaluate(-1,1);  
 }
 
 template<class R>
-multivariate_laurentpoly<Z> compute_khp(const knot_diagram& k, bool reduced) {
-  cube<R> c(kd,reduced);
-  ptr<const module<R> > C = c.khC;
-  mod_map<R> d = c.compute_d(1, 0, 0, 0, 0);
-  chain_complex_simplifier<R> s(C, d, maybe<int>(1), maybe<int>(0));
-  C = s.new_C;
-  d = s.new_d;
-  multivariate_laurentpoly<Z> khp;
-  for(uint i = 1; i <= C->dim(); ++i) {
-    grading gr = C->generator_grading(i);
-    multivariate_laurentpoly<Z> p1 = multivariate_laurentpoly<Z>(1, VARIABLE, 0, gr.h);
-    multivariate_laurentpoly<Z> p2 = multivariate_laurentpoly<Z>(1, VARIABLE, 1, gr.q);
-    multivariate_laurentpoly<Z> p3 = p1 * p2;
-    khp += p3;
+int compute_s_inv(knot_diagram& kd) {
+  unsigned m = kd.num_components ();
+  if (m != 1) {
+    fprintf (stderr, "error: s-invariant only defined for knots\n");
+    exit (EXIT_FAILURE);
   }
-  return khp;
+      
+  cube<R> c (kd, 0);
+  ptr<const module<R> > C = c.khC;
+      
+  mod_map<R> d = c.compute_d (1, 0, 0, 0, 0);
+  for (unsigned i = 1; i <= kd.n_crossings; i ++)
+    d = d + c.H_i (i);
+  assert (d.compose (d) == 0);
+      
+  int k = 0;
+  for (;;) {
+    chain_complex_simplifier<R> s (C, d,
+				   maybe<int> (1),
+				   maybe<int> (2*k));
+    C = s.new_C;
+    d = s.new_d;
+    k ++;
+	  
+    if (d == 0)
+      break;
+  }
+      
+  assert (C->dim () == 2);
+  grading gr1 = C->generator_grading (1),
+    gr2 = C->generator_grading (2);
+  C->free_poincare_polynomial().display_self();
+  int qmin = gr1.q,
+    qmax = gr2.q;
+  if (qmax < qmin)
+    std::swap (qmin, qmax);
+      
+  assert (qmax == qmin + 2);
+  return qmin + 1;
 }
 
+void check_periodicity(std::string out_file) {
+  if(period == 2 || period == 3) {
+    std::cout << "Sorry, the criteria don't work for period "
+	      << period << "...\n";
+    exit(EXIT_FAILURE);
+  }
+  auto result = std::find(primes_list.begin(), primes_list.end(), period);
+  if(result == primes_list.end()) {
+    std::cout << "For now you can only check periodicity for primes up to 31..." << "\n";
+    exit(EXIT_FAILURE);
+  }
+  std::ofstream out(out_file);
+  
+  if(periodicity_test == "all") {
+    std::cout << "I will perform both test..." << "\n";
+  }
+  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;
+      else
+	std::cout << pcc(compute_jones(kd)) << 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);
+    if(out_file.size() != 0)
+      out << pc(period) << std::endl;
+    else
+      std::cout << pc(period) << std::endl;
+  }
+  else {
+    std::cout << "Sorry, I don't recognize this option..." << "\n";
+    exit(EXIT_FAILURE);
+  }
+}
+
+
+// template<class R>
+// void check_periodicity_criterion(knot_diagram kd, int prime = 5, Test_type t){
+  // using monomial = multivariate_laurent_monomial;
+  // using polynomial = multivariate_laurentpoly<Z>;
+  // unsigned m = kd.num_components();
+  // if(m != 1) {
+  //   std::cerr << "Error: for now this only works for knots\n";
+  //   exit(EXIT_FAILURE);
+  // }
+  // polynomial khp, lee_p;
+  // // Adding braces so that the cube can be destroyed
+  // // when its not needed anymore
+  // { 
+  //   cube<R> c(kd,0);
+  //   ptr<const module<R>> C = c.khC;
+  //   mod_map<R> d = c.compute_d(1, 0, 0, 0, 0);
+  //   for(unsigned i = 1; i <= kd.n_crossings; i++)
+  //     d = d + c.H_i(i);
+  //   assert (d.compose(d) == 0);
+
+  //   // computing Khovanov homology
+  //   if(verbose)
+  //     std::cout << "Computing Khovanov polynomial..." << "\n";
+  //   chain_complex_simplifier<R> s(C, d, maybe<int>(1), maybe<int>(0));
+  //   C = s.new_C;
+  //   d = s.new_d;
+  //   khp = C->free_poincare_polynomial();
+
+  //   // computing Lee homology
+  //   if(verbose)
+  //     std::cout << "Computing Lee polynomial..." << "\n";
+  //   chain_complex_simplifier<R> s1(C, d, maybe<int>(1), maybe<int>(2));
+  //   C = s1.new_C;
+  //   d = s1.new_d;
+  //   assert(C->dim() == 2);
+  //   lee_p = C->free_poincare_polynomial();
+  // }
+  // if(verbose) {
+  //   std::cout << "Khovanov polynomial = "
+  // 	      << khp << "\n"
+  // 	      << "Lee polynomial = "
+  // 	      << lee_p << "\n";
+  // }
+  // periodicity_checker pc(khp, lee_p, prime, t);
+  // pc();
+// }
+
 template<class R> void
 compute_invariant ()
 {
@@ -413,48 +570,10 @@ compute_invariant ()
       tex_footer ();
     }
   else if (!strcmp (invariant, "s"))
-    {
-      unsigned m = kd.num_components ();
-      if (m != 1)
-	{
-	  fprintf (stderr, "error: s-invariant only defined for knots\n");
-	  exit (EXIT_FAILURE);
-	}
-      
-      cube<R> c (kd, 0);
-      ptr<const module<R> > C = c.khC;
-      
-      mod_map<R> d = c.compute_d (1, 0, 0, 0, 0);
-      for (unsigned i = 1; i <= kd.n_crossings; i ++)
-	d = d + c.H_i (i);
-      assert (d.compose (d) == 0);
-      
-      int k = 0;
-      for (;;)
-	{
-	  chain_complex_simplifier<R> s (C, d,
-					 maybe<int> (1), maybe<int> (2*k));
-	  C = s.new_C;
-	  d = s.new_d;
-	  k ++;
-	  
-	  if (d == 0)
-	    break;
-	}
-      
-      assert (C->dim () == 2);
-      grading gr1 = C->generator_grading (1),
-	gr2 = C->generator_grading (2);
-      
-      int qmin = gr1.q,
-	qmax = gr2.q;
-      if (qmax < qmin)
-	std::swap (qmin, qmax);
-      
-      assert (qmax == qmin + 2);
-      
-      fprintf (outfp, "s(%s; %s) = %d\n", knot, field, qmin + 1);
-    }
+  {
+    int s = compute_s_inv<R>(kd);  
+    fprintf (outfp, "s(%s; %s) = %d\n", knot, field, s);
+  }
   else 
     {
       fprintf (stderr, "error: unknown invariant %s\n", invariant);
@@ -617,69 +736,73 @@ main (int argc, char **argv)
   
   const char *file = 0;
   
-  for (int i = 1; i < argc; i ++)
-    {
-      if (argv[i][0] == '-')
-	{
-	  if (strcmp (argv[i], "-r") == 0)
-	    reduced = 1;
-	  else if (strcmp (argv[i], "-h") == 0)
-	    {
-	      usage ();
-	      exit (EXIT_SUCCESS);
-	    }
-		else if (strcmp (argv[i], "-demo") == 0)
-    {
-      run_demo();
-      exit (EXIT_SUCCESS);
-    }	  
-	  else if (!strcmp (argv[i], "-v"))
-	    verbose = 1;
-	  else if (!strcmp (argv[i], "-f"))
-	    {
-	      i ++;
-	      if (i == argc)
-		{
-		  fprintf (stderr, "error: missing argument to option `-f'\n");
-		  exit (EXIT_FAILURE);
-		}
-	      field = argv[i];
-	    }
-	  else if (!strcmp (argv[i], "-o"))
-	    {
-	      i ++;
-	      if (i == argc)
-		{
-		  fprintf (stderr, "error: missing argument to option `-o'\n");
-		  exit (EXIT_FAILURE);
-		}
-	      file = argv[i];
-	    }
-	  else
-	    {
-	      fprintf (stderr, "error: unknown argument `%s'\n", argv[1]);
-	      fprintf (stderr, "  use -h for usage\n");
-	      exit (EXIT_FAILURE);
-	    }
+  for (int i = 1; i < argc; i ++) {
+    if (argv[i][0] == '-') {
+      if (strcmp (argv[i], "-r") == 0)
+	reduced = 1;
+      else if (strcmp (argv[i], "-h") == 0) {
+        usage ();
+        exit (EXIT_SUCCESS);
+      }
+      else if (strcmp (argv[i], "-demo") == 0) {
+	run_demo();
+	exit (EXIT_SUCCESS);
+      }	  
+      else if (!strcmp (argv[i], "-v"))
+	verbose = 1;
+      else if (!strcmp (argv[i], "-f")) {
+	i ++;
+	if (i == argc) {
+	  fprintf (stderr, "error: missing argument to option `-f'\n");
+	  exit (EXIT_FAILURE);
 	}
-      else
-	{
-	  if (knot)
-	    {
-	      fprintf (stderr, "error: too many arguments\n");
-	      fprintf (stderr, "  use -h for usage\n");
-	      exit (EXIT_FAILURE);
-	    }
-	  else if (invariant)
-	    knot = argv[i];
-	  else
-	    {
-	      assert (invariant == 0);
-	      invariant = argv[i];
-	    }
+	field = argv[i];
+      }
+      else if (!strcmp (argv[i], "-o")) {
+	i ++;
+	if (i == argc) {
+	  fprintf (stderr, "error: missing argument to option `-o'\n");
+	  exit (EXIT_FAILURE);
 	}
+	file = argv[i];
+      }
+      else if(!strcmp (argv[i], "-p")) {
+	i++;
+	if(i == argc) {
+	  fprintf (stderr, "error: missing argument to option `-o'\n");
+	  exit (EXIT_FAILURE);
+	}
+	period = std::stoi(argv[i]);
+      }
+      else if (!strcmp(argv[i], "-t")) {
+	i++;
+	if(i == argc) {
+	  fprintf (stderr, "error: missing argument to option `-o'\n");
+	  exit (EXIT_FAILURE);
+	}
+	periodicity_test = argv[i];
+      }
+      else {
+	fprintf (stderr, "error: unknown argument `%s'\n", argv[1]);
+	fprintf (stderr, "  use -h for usage\n");
+	exit (EXIT_FAILURE);
+      }
     }
-  
+    else {
+      if (knot) {
+	fprintf (stderr, "error: too many arguments\n");
+	fprintf (stderr, "  use -h for usage\n");
+	exit (EXIT_FAILURE);
+      }
+      else if (invariant)
+	knot = argv[i];
+      else {
+	assert (invariant == 0);
+	invariant = argv[i];
+      }
+    }
+  }
+ 
   if (!knot)
     {
       fprintf (stderr, "error: too few arguments, <invariant> or <knot> missing\n");
@@ -687,15 +810,13 @@ main (int argc, char **argv)
       exit (EXIT_FAILURE);
     }
   
-  if (file)
-    {
-      outfp = fopen (file, "w");
-      if (!outfp)
-	{
-	  stderror ("fopen: %s", file);
-	  exit (EXIT_FAILURE);
-	}
+  if (file) {
+    outfp = fopen (file, "w");
+    if (!outfp) {
+      stderror ("fopen: %s", file);
+      exit (EXIT_FAILURE);
     }
+  }
   else
     outfp = stdout;
   
@@ -740,36 +861,31 @@ main (int argc, char **argv)
       compute_gss ();
     }
   else if(!strcmp(invariant, "jones")) {
-    multivariate_laurentpoly<Z> jones_pol = compute_jones(kd, reduced);
-    printf("Jones polynomial of %s is equal to (coefficients in %s):\n", knot, field);
-    jones_pol.display_self();
+    std::cout << "Jones polynomial of " << knot << " = " << compute_jones(kd, reduced) << "\n";
   }
-  else if(!strcmp(invariant, "przytycki_cong")) {
-    // need to add a commandline switch to be able to read
-    // symmetry order
-    multivariate_laurentpoly<Z> jones_pol = compute_jones(kd, reduced);
-    multivariate_laurentpoly<Z> inv_jones_pol = invert_variable(jones_pol,1);
-    jones_pol.display_self();
-    inv_jones_pol.display_self();
-    multivariate_laurentpoly<Z> diff = jones_pol - inv_jones_pol;
-    diff.display_self();
-    check_przytycki_cong(diff, 3);
+  else if(!strcmp(invariant, "periodicity")) {
+    check_periodicity((file ? std::string(file) : std::string()));
   }
   else if(!strcmp(invariant, "khp")) {
     multivariate_laurentpoly<Z> khp;
     if(!strcmp(field, "Z2"))
       khp = compute_khp<Z2>(kd, reduced);
     else if(!strcmp(field, "Z3"))
-      khp = compute_khp<Zp<3> >(kd, reduced);
+      khp = compute_khp<Zp<3>>(kd, reduced);
+    else if(!strcmp(field, "Z5"))
+      khp = compute_khp<Zp<5>>(kd, reduced);
+    else if (!strcmp(field, "Z7"))
+      khp = compute_khp<Zp<7>>(kd,reduced);
     else if(!strcmp(field, "Q"))
       khp = compute_khp<Q>(kd, reduced);
     else
     {
-      fprintf (stderr, "error: unknown field %s\n", field);
+      std::cerr << "Unknown field: " << field << std::endl;
       exit (EXIT_FAILURE);
     }
-    printf("Khovanov polynomial of %s is equal to (coefficients %s):\n", knot, field);
-    khp.display_self();
+    std::cout << "Khovanov polynomial (coefficients in " << field
+	      << ") of " << knot <<  " = " << std::endl
+	      << khp << std::endl;
   }
   else
     {
diff --git a/knot_diagram.h b/knot_diagram.h
index 8cd3808..e05719a 100644
--- a/knot_diagram.h
+++ b/knot_diagram.h
@@ -1,3 +1,7 @@
+#ifndef _KNOTKIT_KNOT_DIAGRAM_H
+#define _KNOTKIT_KNOT_DIAGRAM_H
+
+#include <planar_diagram.h>
 
 // for building knot_diagram
 inline unsigned edge_from_ept (unsigned e)
@@ -183,3 +187,5 @@ class knot_diagram
   void show_self () const;
   void display_self () const;
 };
+
+#endif // _KNOTKIT_KNOT_DIAGRAM_H
diff --git a/knotkit.h b/knotkit.h
index b1974e6..27f4a39 100644
--- a/knotkit.h
+++ b/knotkit.h
@@ -1,3 +1,5 @@
+#ifndef _KNOTKIT_KNOTKIT_H
+#define _KNOTKIT_KNOTKIT_H
 
 // includes lib.h
 #include <algebra/algebra.h>
@@ -123,3 +125,5 @@ resolution_diagram parse_resolution_diagram (const char *s);
 
 // 11 <= n <= 15
 basedvector<basedvector<unsigned, 1>, 1> mutant_knot_groups (unsigned n);
+
+#endif // _KNOTKIT_KNOTKIT_H
diff --git a/lib/bitset.h b/lib/bitset.h
index 48cb575..f690309 100644
--- a/lib/bitset.h
+++ b/lib/bitset.h
@@ -1,3 +1,8 @@
+#ifndef _KNOTKIT_LIB_BITSET_H
+#define _KNOTKIT_LIB_BITSET_H
+
+#include <lib/lib.h>
+#include <lib/vector.h>
 
 typedef uint64 word_t;
   
@@ -136,3 +141,5 @@ bitset::operator % (unsigned i) const
   unsigned b = (i - 1) & word_bit_mask;
   return word_bittest (v[w], b + 1);
 }
+
+#endif // _KNOTKIT_LIB_BITSET_H
diff --git a/lib/lib.h b/lib/lib.h
index 4537ade..2e0c1fb 100644
--- a/lib/lib.h
+++ b/lib/lib.h
@@ -1,3 +1,5 @@
+#ifndef _KNOTKIT_LIB_LIB_H
+#define _KNOTKIT_LIB_LIB_H
 
 #include <assert.h>
 #include <stdlib.h>
@@ -235,3 +237,5 @@ public:
   map<unsigned, unsigned> io_id_id;
   
 };
+
+#endif // _KNOTKIT_LIB_LIB_H
diff --git a/lib/vector.h b/lib/vector.h
index 06f4e52..7650733 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -1,3 +1,7 @@
+#ifndef _KNOTKIT_LIB_VECTOR_H
+#define _KNOTKIT_LIB_VECTOR_H
+
+#include <lib/lib.h>
 
 template<class T>
 class vector
@@ -386,3 +390,5 @@ public:
   unsigned lower_bound (const T &v) const { return vector<T>::lower_bound (v) + B; }
   unsigned upper_bound (const T &v) const { return vector<T>::upper_bound (v) + B; }
 };
+
+#endif // _KNOTKIT_LIB_VECTOR_H
diff --git a/periodicity.cpp b/periodicity.cpp
new file mode 100644
index 0000000..6d98527
--- /dev/null
+++ b/periodicity.cpp
@@ -0,0 +1,158 @@
+#include <periodicity.h>
+#include <simplify_chain_complex.h>
+
+void Kh_periodicity_checker::compute_knot_polynomials(knot_diagram& kd) {
+
+  unsigned m = kd.num_components ();
+  if (m != 1) {
+    std::cerr << "warning: this implementation of the criterion works for knots only...";
+    exit (EXIT_FAILURE);
+  }
+      
+  cube<Z2> c (kd, 0);
+  ptr<const module<Z2> > C = c.khC;
+      
+  mod_map<Z2> d = c.compute_d (1, 0, 0, 0, 0);
+  for (unsigned i = 1; i <= kd.n_crossings; i ++)
+    d = d + c.H_i (i);
+  assert (d.compose (d) == 0);
+
+  // computing Khovanov homology
+  if(verbose)
+    std::cout << "Computing Khovanov homology" << std::endl;
+  {
+    chain_complex_simplifier<Z2> s (C, d, maybe<int>(1), maybe<int>(0));
+    C = s.new_C;
+    d = s.new_d;
+    khp = C->free_poincare_polynomial();
+    if(verbose)
+      std::cout << "KhP = " << khp << "\n";
+  }
+  
+  // computing Lee homolgy
+  if(verbose)
+    std::cout << "Computing Lee homology" << std::endl;
+  {
+    chain_complex_simplifier<Z2> s(C, d, maybe<int>(1), maybe<int>(2));
+    C = s.new_C;
+    d = s.new_d;
+    leep = C->free_poincare_polynomial();
+    if(d != 0) {
+      std::cout << "For now, you can only use this criterion on Kh-thin knots." << std::endl;
+      exit(EXIT_FAILURE);
+    }
+    if(verbose) {
+      std::cout << "LeeP = " << leep << "\n";
+    }
+  }
+}
+
+void Kh_periodicity_checker::compute_quot() {
+  polynomial diff = khp - leep;
+  while(diff != 0) {
+    pair<monomial, Z> m = diff.head();
+    if(m.first.m[1] == 1) {
+      pair<monomial, Z> m1 = diff.tail();
+      while(m1.first.m.card() == 1 && m1.first.m[2]) {
+	quot += polynomial(m1.second, m1.first);
+	polynomial p = polynomial(m1.second, m1.first) * mul;
+	diff -= p;
+	m1 = diff.tail();
+      }
+    }
+    quot += polynomial(m.second, m.first);
+    polynomial p = polynomial(m.second, m.first) * mul;
+    diff -= p;
+  }
+}
+
+std::pair<multivariate_laurentpoly<Z>, multivariate_laurentpoly<Z>>
+Kh_periodicity_checker::compute_quotient_and_remainder(const polynomial& quot,
+						       int period) const {
+  polynomial quotient, remainder;
+  for(map<monomial, Z>::const_iter i = quot.coeffs; i; i++) {
+    std::tuple<Z,Z> div = i.val().divide_with_remainder(period - 1);
+    quotient += polynomial(std::get<0>(div), i.key());
+    remainder += polynomial(std::get<1>(div), i.key());
+  }
+  if(verbose) {
+    std::cout << "Decomposition of Khp = " << std::endl
+	      << leep << " + ("
+	      << mul << ") * ("
+	      << remainder;
+    if(quotient != 0) {
+      std::cout << " + " << (period - 1)
+		<< " * (" << quotient
+		<< ")";
+    }
+    std::cout << ")" << std::endl;
+   
+  }
+  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 {
+  periodic_congruence_checker<Z> pcc(period);
+  polynomial t = leep + mul * r;
+  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;
+    polynomial m = mul;
+    return any_of(candidates.begin(), candidates.end(),
+		  [&pcc, &t, &m, &period](const polynomial& p){
+		    return pcc((t + Z(period - 1) * m * p).evaluate(-1,1));
+      });
+  }
+}
+
+std::string Kh_periodicity_checker::operator () (int period) const {
+  std::ostringstream out;
+  std::pair<polynomial, polynomial> q_r = compute_quotient_and_remainder(quot, period);
+  bool res = check(std::get<0>(q_r), std::get<1>(q_r), period);
+  out << knot << ": period = " << period << ": "
+      << (res ? "Maybe" : "No");
+  return out.str();
+}
diff --git a/periodicity.h b/periodicity.h
new file mode 100644
index 0000000..616b5bb
--- /dev/null
+++ b/periodicity.h
@@ -0,0 +1,124 @@
+#ifndef _KNOTKIT_PERIODICITY_H
+#define _KNOTKIT_PERIODICITY_H
+
+#include <knotkit.h>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <utility>
+#include <list>
+#include <algorithm>
+
+extern bool verbose;
+extern const char* knot;
+
+extern std::string periodicity_test;
+
+const std::vector<int> primes_list = {5, 7, 11, 13, 17, 19};
+
+enum class Test_type {
+  Przytycki_criterion,
+  BKP_criterion,
+  all
+};
+
+template<class T>
+class periodic_congruence_checker {
+public:
+  using polynomial = multivariate_laurentpoly<T>;
+  using monomial = multivariate_laurent_monomial;
+  
+  int prime;
+  unsigned index;
+
+  polynomial prepare_polynomial(const polynomial& pol) const {
+    return (pol - invert_variable(pol, index));
+  }
+
+  bool reduce(const polynomial& pol) const;
+
+public:
+  periodic_congruence_checker(int pprime = 5,
+			      unsigned ind = 2) :
+    prime(pprime),
+    index(ind)
+    {}
+
+  ~periodic_congruence_checker() {};
+
+  bool operator() (const polynomial& pol) const {
+    return reduce(prepare_polynomial(pol));
+  }
+};
+
+template<class T> 
+bool periodic_congruence_checker<T>::reduce(const multivariate_laurentpoly<T>& pol) const {
+  polynomial res;
+  for(typename map<monomial, T>::const_iter i = pol.coeffs; i; i++) {
+    int c = i.key().m[index] % (2 * prime);
+    if(c < 0)
+      c += (2 * prime);
+    monomial mon = monomial(VARIABLE, index, c);
+    res += polynomial(i.val(), mon);
+  }
+  if(verbose)
+    std::cout << "reduced = " << res << "\n";
+  return res == 0;
+}
+
+template<unsigned p>
+class Przytycki_periodicity_checker {
+  using polynomial = multivariate_laurentpoly<Zp<p>>;
+  using monomial = multivariate_laurent_monomial;
+
+  periodic_congruence_checker<Zp<p>> cong_checker;
+
+ public:
+  Przytycki_periodicity_checker() : cong_checker(p) {}
+    
+  ~Przytycki_periodicity_checker() {}
+  
+  std::string operator() (const polynomial& pol) const {
+    std::ostringstream out;
+    out << knot << ": period = " << p
+	<< ": "
+	<< (cong_checker(pol) ? "Maybe" : "No");
+    return out.str();
+  }
+};
+
+class Kh_periodicity_checker {
+  using polynomial = multivariate_laurentpoly<Z>;
+  using monomial = multivariate_laurent_monomial;
+
+  unsigned ev_index;
+  unsigned index;
+
+  polynomial khp, leep, quot;
+  polynomial mul;
+
+  void compute_knot_polynomials(knot_diagram& kd);
+  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;
+  bool check(const polynomial& q, const polynomial& r, int period) const;
+
+ public:
+  Kh_periodicity_checker(knot_diagram& kd) {
+    ev_index = 1;
+    index = 2;
+    mul = polynomial(Z(1))
+      + polynomial(1, VARIABLE, ev_index) *
+      polynomial(1, VARIABLE, index, 2);
+
+    compute_knot_polynomials(kd);
+    compute_quot();
+  }
+
+  ~Kh_periodicity_checker() {}
+
+  std::string operator () (int period) const;
+};
+
+#endif // _KNOTKIT_PERIODICITY_H
diff --git a/planar_diagram.h b/planar_diagram.h
index 6c32c3a..e39a721 100644
--- a/planar_diagram.h
+++ b/planar_diagram.h
@@ -1,7 +1,11 @@
-
+#ifndef _KNOTKIT_PLANAR_DIAGRAM_H
+#define _KNOTKIT_PLANAR_DIAGRAM_H
 /* Planar diagram of a knot.  For details, see:
    http://katlas.org/wiki/Planar_Diagrams */
 
+#include <lib/vector.h>
+#include <string>
+
 class planar_diagram
 {
 public:
@@ -27,3 +31,5 @@ public:
   void show_self () const { printf ("planar_diagram %s", name.c_str ()); }
   void display_self () const;
 };
+
+#endif // _KNOTKIT_PLANAR_DIAGRAM_H
diff --git a/simplify_chain_complex.h b/simplify_chain_complex.h
index 189e48d..a302532 100644
--- a/simplify_chain_complex.h
+++ b/simplify_chain_complex.h
@@ -1,3 +1,5 @@
+#ifndef _KNOTKIT_SIMPLIFY_CHAIN_COMPLEX_H
+#define _KNOTKIT_SIMPLIFY_CHAIN_COMPLEX_H
 
 template<class R> class simplified_complex_generators
 {
@@ -247,3 +249,5 @@ chain_complex_simplifier<R>::chain_complex_simplifier (ptr<const module<R> > C_,
   assert (new_d.compose (pi) == pi.compose (d));
   assert (pi.compose (iota) == mod_map<R> (new_C, 1));
 }
+
+#endif // _KNOTKIT_SIMPLIFY_CHAIN_COMPLEX_H