diff --git a/Makefile b/Makefile
index d79926f..c154cff 100644
--- a/Makefile
+++ b/Makefile
@@ -12,7 +12,7 @@ INCLUDES = -I. -I/opt/local/include
 
 # OPTFLAGS = -g
 OPTFLAGS = -O2 -g
-# OPTFLAGS = -O2 -g -DNDEBUG
+# OPTFLAGS = -O2 -DNDEBUG
 
 LDFLAGS = -L/opt/local/lib
 
diff --git a/algebra/multivariate_laurentpoly.h b/algebra/multivariate_laurentpoly.h
index 6c9d548..8b1e04a 100644
--- a/algebra/multivariate_laurentpoly.h
+++ b/algebra/multivariate_laurentpoly.h
@@ -365,7 +365,7 @@ class multivariate_laurentpoly
   void check () const;
 #endif
 
-  multivariate_laurentpoly evaluate(T val, unsigned index) const;
+  multivariate_laurentpoly& evaluate(T val, unsigned index);
   
   static void show_ring ()
   {
@@ -515,15 +515,6 @@ std::ostream& operator << (std::ostream& os, const multivariate_laurentpoly<T>&
   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
 
 // function below inverts every occurence of a variable x_index
@@ -550,22 +541,34 @@ multivariate_laurentpoly<T> invert_variable(const multivariate_laurentpoly<T>& p
 
 template<class T>
 multivariate_laurentpoly<T>
-multivariate_laurentpoly<T>::evaluate(T val, unsigned index) const {
+evaluate_with_copy(multivariate_laurentpoly<T> pol,
+		   T val, unsigned index) {
   using polynomial = multivariate_laurentpoly<T>;
   using monomial = multivariate_laurent_monomial;
   polynomial res;
-  for(typename map<monomial, T>::const_iter i = coeffs; i; i++) {
+  for(typename map<monomial, T>::const_iter i = pol.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;
+      monomial mon;
+      for(map<unsigned, int>::const_iter j = i.key().m; j; ++j) {
+	if(j.key() != index)
+	  mon *= monomial(VARIABLE, j.key(), j.val());
+      }
+      res += polynomial(i.val() * pow(val, exp), mon);
+    }
+    else {
+      monomial mon(COPY, i.key());
+      res += polynomial(i.val(), mon);
     }
-    else
-      res += polynomial(i.val(), i.key());
   }
   return res;
 }
 
+template<class T>
+multivariate_laurentpoly<T>&
+multivariate_laurentpoly<T>::evaluate(T val, unsigned index) {
+  *this = evaluate_with_copy<T>(*this, val, index);
+  return *this;
+}
+
 #endif // _KNOTKIT_ALGEBRA_MULTIVARIATE_LAURENPOLY_H
diff --git a/kk.cpp b/kk.cpp
index 441f3f9..01242dd 100644
--- a/kk.cpp
+++ b/kk.cpp
@@ -3,6 +3,7 @@
 #include <fstream>
 #include <vector>
 #include <utility>
+#include <cctype>
 
 const char *program_name;
 
@@ -283,7 +284,7 @@ multivariate_laurentpoly<Z> compute_khp(const knot_diagram& k, bool reduced = fa
 }
 
 multivariate_laurentpoly<Z> compute_jones(const knot_diagram& k, bool reduced = false) {
-  return compute_khp<Z2>(k, reduced).evaluate(-1,1);  
+  return compute_khp<Z2>(k, reduced).evaluate(-1, 1);  
 }
 
 template<class R>
@@ -330,13 +331,63 @@ 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));
+    Kh_periodicity_checker Kh_pc(kd, std::string(knot));
     for(auto& p : primes_list) {
-      std::cout << "Przytycki's criterion: "
-		<< P_pc(p) << std::endl
-		<< "Kh criterion: "
-		<< Kh_pc(p) << std::endl;
+      std::cout << "Kh criterion: "
+    		<< Kh_pc(p) << std::endl;
+    }
+  }
+  else if(periodicity_test == "all_seq") {
+    std::string k(knot);
+    // first check whether the number of crossings is bigger or lower than 10
+    if(isdigit(k[1])) {
+      // at least 10 crossings
+      if(k[1] == '0') {
+	// ten crossings
+	int num_cr = 10;
+	int knot_index = stoi(k.substr(3));
+	for(int i = knot_index; i < rolfsen_crossing_knots(num_cr); i++) {
+	  std::string knot_name = std::to_string(num_cr) + "_" + std::to_string(i);
+	  knot_diagram kd_temp = parse_knot(knot_name.c_str());
+	  kd.marked_edge = 1;
+	  Kh_periodicity_checker Kh_pc(kd_temp, knot_name);
+	  for(auto& p : primes_list) {
+	    std::cout << "Kh criterion: "
+		      << Kh_pc(p) << std::endl;
+	  }
+	}
+      }
+      else {
+	int num_cr = stoi(k.substr(0,2));
+	int knot_index = stoi(k.substr(3));
+	char alt = k[2];
+	bool alternating = (alt == 'a' ? true : false);
+	for(int i = knot_index; i <= htw_knots(num_cr, alternating); i++) {
+	  std::string knot_name = std::to_string(num_cr) + alt + std::to_string(i);
+	  knot_diagram kd_temp = parse_knot(knot_name.c_str());
+	  kd.marked_edge = 1;
+	  Kh_periodicity_checker Kh_pc(kd_temp, knot_name);
+	  for(auto& p : primes_list) {
+	    std::cout << "Kh criterion: "
+		      << Kh_pc(p) << std::endl;
+	  }
+	}
+      }
+    }
+    else {
+      // at most nine crossings
+      int num_cr = stoi(k.substr(0, 1));
+      int knot_index = stoi(k.substr(2));
+      for(int i = knot_index; i <= rolfsen_crossing_knots(num_cr); i++) {
+	std::string knot_name = std::to_string(num_cr) + "_" + std::to_string(i);
+	knot_diagram kd_temp = parse_knot(knot_name.c_str());
+	kd.marked_edge = 1;
+	Kh_periodicity_checker Kh_pc(kd_temp, knot_name);
+	for(auto& p : primes_list) {
+	  std::cout << "Kh criterion: "
+		    << Kh_pc(p) << std::endl;
+	}
+      }
     }
   }
   else {
@@ -360,7 +411,7 @@ void check_periodicity(std::string out_file) {
       std::cout << P_pc(period) << std::endl;
   }
   else if(periodicity_test == "Kh") {
-    Kh_periodicity_checker Kh_pc(kd);
+    Kh_periodicity_checker Kh_pc(kd, std::string(knot));
     if(out_file.size() != 0)
       out << Kh_pc(period) << std::endl;
     else
diff --git a/periodicity.cpp b/periodicity.cpp
index 085324c..468a523 100644
--- a/periodicity.cpp
+++ b/periodicity.cpp
@@ -28,6 +28,7 @@ bool Przytycki_periodicity_checker::check(int period) const {
     return pcc(jones_pol);
   }
   }
+  return false;
 }
 
 std::string Przytycki_periodicity_checker::operator () (int period) const {
@@ -37,114 +38,6 @@ std::string Przytycki_periodicity_checker::operator () (int period) const {
   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 ();
@@ -241,49 +134,123 @@ Kh_periodicity_checker::compute_quotient_and_remainder(const polynomial& quot,
   return std::make_pair(quotient, remainder);
 }
 
+std::map<multivariate_laurentpoly<Z>, std::pair<Z,Z>>
+Kh_periodicity_checker::compute_bounds(const polynomial& p, int period) const {
+  std::map<polynomial, std::pair<Z, Z>> bounds;
+  periodic_congruence_checker<Z> pcc(period);
+  for(map<monomial, Z>::const_iter i = p.coeffs; i; ++i) {
+    monomial mon;
+    int exp = 0;
+    if(i.key().m % ev_index) {
+      exp = i.key().m[ev_index];
+      for(map<unsigned, int>::const_iter j = i.key().m; j; ++j) {
+	if(j.key() != ev_index) {
+	  int v = j.val() % (2 * period);
+	  if(v < 0) v += (2 * period);
+	  mon *= monomial(VARIABLE, j.key(), v);
+	}
+      }
+    }
+    else {
+      for(map<unsigned, int>::const_iter j = i.key().m; j; ++j) {
+	int v = j.val() % (2 + period);
+	if (v < 0) v += (2 * period);
+	mon *= monomial(VARIABLE, j.key(), v);
+      }
+    }
+    // std::cout << polynomial(i.val() * pow(-1, exp), mon) << "\n";
+    Z v_temp = i.val() * pow(-1, exp);
+    polynomial p_temp = (polynomial(1, mon) * mul).evaluate(-1, ev_index);
+    p_temp = pcc.reduce(p_temp - invert_variable(p_temp, index));
+    if(v_temp >= 0)
+      bounds[p_temp].second += (v_temp * period);
+    else
+      bounds[p_temp].first += (v_temp * period);
+  }
+  
+  // for(std::map<polynomial, std::pair<Z,Z>>::iterator mi = bounds.begin(); mi != bounds.end(); ++mi) {
+  //   std::cout << "Monomial: " << mi->first << "\n";
+  //   std::cout << "Max: " << std::get<1>(mi->second)
+  // 	      << ", Min: " << std::get<0>(mi->second) << "\n";
+  // }
+  return bounds;
+}
+
+std::vector<multivariate_laurentpoly<Z>>
+Kh_periodicity_checker::compute_basis_polynomials(int period) const {
+  std::vector<polynomial> res;
+  periodic_congruence_checker<Z> pcc(period);
+  for(int i = 1; i < period; i += 2) {
+    res.push_back(pcc.reduce(get_basis_polynomial(i)));
+  }
+  return res;
+}
+
+multivariate_laurentpoly<Z> Kh_periodicity_checker::get_basis_polynomial(monomial mon) const {
+  return (polynomial(Z(1), mon) * mul).evaluate(-1, ev_index) -
+    invert_variable((polynomial(Z(1), mon) * mul).evaluate(-1, ev_index), index);
+}
+
 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));
+  polynomial t = (leep + mul * (r - q)).evaluate(-1, ev_index);
+  t = pcc.reduce(t - invert_variable(t, index));
+  if(pcc(t)) {
+    return true;
   }
-  if(verbose)
-    std::cout << "Checking congruences..."; 
-  polynomial_iterator<Z> pi(q);
-  polynomial_iterator<Z> pi_end(q, polynomial_iterator<Z>::start_pos::end);
-  Z count = pi.get_count();
-  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;
+  else if(q == 0)
+    return false;
+  // std::cout << t << std::endl;
+  // std::cout << q << "\n";
+  std::map<polynomial, std::pair<Z,Z>> bounds = compute_bounds(q, period);
+  for(std::map<polynomial, std::pair<Z,Z>>::iterator it = bounds.begin();
+      it != bounds.end(); ++it) {
+    polynomial mon = it->first;
+  }
+  std::vector<polynomial> basis_polynomials = compute_basis_polynomials(period);
+  polynomial p = pcc.reduce(get_basis_polynomial(2 * period - 1));
+  for(Z i = bounds[p].first; i <= bounds[p].second; i += 5) {
+    polynomial p_temp = t + polynomial(i, VARIABLE, index, 0) * p;
+    // std::cout << "i = " << i << "\n";
+    // std::cout << "p_temp = " << p_temp << "\n";
+    if(p_temp == 0)
       return true;
+    for(std::vector<polynomial>::iterator it = basis_polynomials.begin(); it != basis_polynomials.end(); ++it) {
+      pair<monomial, Z> m = p_temp.coeffs.head();
+      monomial mon = m.first;
+      Z c = m.second;
+      polynomial pp = pcc.reduce(get_basis_polynomial(mon));
+      if(pp == *it) {
+	if(c < bounds[pp].first || c > bounds[pp].second)
+	  break;
+	else {
+	  // std::cout << "pp = " << pp << "\n";
+	  p_temp -= polynomial(c, VARIABLE, index, 0) * pp;
+	  // std::cout << "p_temp = " << p_temp << "\n";
+	  if(p_temp == 0)
+	    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 {
   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");
+  // first check Przytycki's criterion
+  Przytycki_periodicity_checker P_pc(evaluate_with_copy<Z>(khp, -1, ev_index));
+  if(!P_pc.check(period)) {
+    out << knot_name << ": period =  " << period << ": No (Przytycki's criterion).";
+  }
+  else {
+    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_name << ": period = " << period << ": "
+	<< (res ? "Maybe" : "No");
+  }
   return out.str();
 }
diff --git a/periodicity.h b/periodicity.h
index 6c83e57..9e37fb7 100644
--- a/periodicity.h
+++ b/periodicity.h
@@ -5,6 +5,7 @@
 #include <sstream>
 #include <string>
 #include <vector>
+#include <utility>
 
 extern bool verbose;
 extern const char* knot;
@@ -30,8 +31,6 @@ protected:
     return pol - inv;
   }
 
-  bool reduce(const polynomial& pol) const;
-
 public:
   periodic_congruence_checker(int pprime = 5,
 			      unsigned ind = invert_index) :
@@ -41,13 +40,17 @@ public:
 
   virtual ~periodic_congruence_checker() {};
 
+  const polynomial reduce(const polynomial& pol) const;
+  
   bool operator() (const polynomial& pol) const {
-    return reduce(prepare_polynomial(pol));
+    return reduce(prepare_polynomial(pol)) == 0;
+    // return reduce(prepare_polynomial(pol)) == 0;
   }
 };
 
 template<class T> 
-bool periodic_congruence_checker<T>::reduce(const multivariate_laurentpoly<T>& pol) const {
+const multivariate_laurentpoly<T>
+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);
@@ -56,9 +59,7 @@ bool periodic_congruence_checker<T>::reduce(const multivariate_laurentpoly<T>& p
     monomial mon = monomial(VARIABLE, index, c);
     res += polynomial(i.val(), mon);
   }
-  // if(verbose)
-  //   std::cout << "reduced = " << res << "\n";
-  return res == 0;
+  return res;
 }
 
 class Przytycki_periodicity_checker {
@@ -67,79 +68,16 @@ class Przytycki_periodicity_checker {
 
   polynomial jones_pol;
 
-  bool check(int period) const;
-
  public:
   Przytycki_periodicity_checker(polynomial j) : jones_pol(j) {}
     
   ~Przytycki_periodicity_checker() {}
+
+  bool check(int period) const;
   
   std::string operator() (int period) const;
 };
 
-template<class T>
-class polynomial_iterator {
-  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;
@@ -150,15 +88,25 @@ class Kh_periodicity_checker {
   polynomial khp, leep, quot;
   polynomial mul;
 
+  std::string knot_name;
+
   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;
+  std::map<polynomial, std::pair<Z,Z>>
+  compute_bounds(const polynomial& p, int period) const;
+  polynomial get_basis_polynomial(int exp) const {
+    return (polynomial(1, VARIABLE, index, exp) * mul).evaluate(-1, ev_index) -
+      invert_variable((polynomial(1, VARIABLE, index, exp) * mul).evaluate(-1, ev_index), index);
+  }
+  polynomial get_basis_polynomial(monomial mon) const;
+  std::vector<polynomial> compute_basis_polynomials(int period) const;
   bool check(const polynomial& q, const polynomial& r, int period) const;
 
  public:
-  Kh_periodicity_checker(knot_diagram& kd) {
+  Kh_periodicity_checker(knot_diagram& kd, std::string knot_n) :
+    knot_name(knot_n) {
     ev_index = 1;
     index = 2;
     mul = polynomial(Z(1))