updated Google hash map support

1 files changed, 711 insertions, 1166 deletions

diff --git a/ripser.cpp b/ripser.cpp
index 97a2158..d2da6e0 100644
--- a/ripser.cpp
+++ b/ripser.cpp
@@ -1,25 +1,28 @@
-#include <vector>
 #include <iostream>
-#include <iomanip>  
 #include <fstream>
-#include <iterator>
-#include <string>
 #include <cassert>
 #include <queue>
+#include <cmath>
+
+#ifdef USE_GOOGLE_HASHMAP
+
 #include <sparsehash/sparse_hash_map>
-#include "prettyprint.hpp"
+template <class Key, class T> class hash_map : public google::sparse_hash_map<Key, T> {
+    public: inline void reserve(size_t hint) { this->resize(hint); } };
+
+#else
+
+#include <unordered_map>
+template <class Key, class T> class hash_map : public std::unordered_map<Key, T> {};
+
+#endif
 
 typedef float value_t;
-//typedef uint16_t value_t;
+// typedef uint16_t value_t;
 
 typedef long index_t;
 typedef long coefficient_t;
 
-//#define PRECOMPUTE_DIAMETERS
-//#define PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION
-
-//#define STORE_DIAMETERS
-
 #define USE_BINARY_SEARCH
 //#define USE_EXPONENTIAL_SEARCH
 
@@ -34,1279 +37,821 @@ typedef long coefficient_t;
 //#define FILE_FORMAT_LOWER_TRIANGULAR_CSV
 
 class binomial_coeff_table {
-    std::vector<std::vector<index_t> > B;
-    index_t n_max, k_max;
-    
+	std::vector<std::vector<index_t>> B;
+	index_t n_max, k_max;
+
 public:
-    binomial_coeff_table(index_t n, index_t k) {
-        n_max = n;
-        k_max = k;
-    
-        B.resize(n + 1);
-        for( index_t i = 0; i <= n; i++ ) {
-            B[i].resize(k + 1);
-            for ( index_t j = 0; j <= std::min(i, k); j++ )
-            {
-                if (j == 0 || j == i)
-                    B[i][j] = 1;
-                else
-                    B[i][j] = B[i-1][j-1] + B[i-1][j];
-            }
-        }
-    }
-    
-    inline index_t operator()(index_t n, index_t k) const {
-        assert(n <= n_max);
-        assert(k <= k_max);
-        return B[n][k];
-    }
+	binomial_coeff_table(index_t n, index_t k) {
+		n_max = n;
+		k_max = k;
+
+		B.resize(n + 1);
+		for (index_t i = 0; i <= n; i++) {
+			B[i].resize(k + 1);
+			for (index_t j = 0; j <= std::min(i, k); j++) {
+				if (j == 0 || j == i)
+					B[i][j] = 1;
+				else
+					B[i][j] = B[i - 1][j - 1] + B[i - 1][j];
+			}
+		}
+	}
+
+	inline index_t operator()(index_t n, index_t k) const {
+		assert(n <= n_max);
+		assert(k <= k_max);
+		return B[n][k];
+	}
 };
 
-//
-// https://comeoncodeon.wordpress.com/2011/10/09/modular-multiplicative-inverse/
-//
-
-// a * (m / a) + m % a = m
-// m % a = -a * (m / a)   (mod m)
-//Dividing by (a * (m % a)):
-// inverse(a) = - (m / a) * inverse(m % a)   (mod m)
-
-std::vector<coefficient_t> multiplicative_inverse_vector (const coefficient_t m) {
-    std::vector<coefficient_t> mod_inverse(m);
-    mod_inverse[1] = 1;
-    for (coefficient_t i = 2; i < m; ++i) {
-        mod_inverse[i] = (-(m/i) * mod_inverse[m % i]) % m + m;
-    }
-    return mod_inverse;
+std::vector<coefficient_t> multiplicative_inverse_vector(const coefficient_t m) {
+	std::vector<coefficient_t> inverse(m);
+	inverse[1] = 1;
+	for (coefficient_t a = 2; a < m; ++a) {
+		// m = a * (m / a) + m % a
+		// Multipying with inverse(a) * inverse(m % a):
+		// 0 = (m / a) * inverse(m % a) + inverse(a) = 0   (mod m)
+		inverse[a] = m - ((m / a) * inverse[m % a]) % m;
+	}
+	return inverse;
 }
 
-template<typename OutputIterator>
-inline OutputIterator get_simplex_vertices( index_t idx, const index_t dim, index_t n, const binomial_coeff_table& binomial_coeff, OutputIterator out  )
-{
-    --n;
-    
-    for( index_t k = dim + 1; k > 0; --k ) {
- 
-        #ifdef USE_BINARY_SEARCH
-        if ( binomial_coeff( n , k ) > idx ) {
-            index_t count;
-            
-            #ifdef USE_EXPONENTIAL_SEARCH
-            for (count = 1; (binomial_coeff( n - count , k ) > idx); count = std::min(count << 1, n));
-            #else
-            count = n;
-            #endif
-            
-            while (count > 0) {
-                index_t i = n;
-                index_t step = count >> 1;
-                i -= step;
-                if (binomial_coeff( i , k ) > idx) {
-                    n = --i;
-                    count -= step + 1;
-                } else count = step;
-            }
-        }
-        #else
-        while( binomial_coeff( n , k ) > idx )
-            --n;
-        #endif
-        
-        assert( binomial_coeff( n , k ) <= idx );
-        assert( binomial_coeff( n + 1, k ) > idx );
-        
-        *out++ = n;
-        
-        idx -= binomial_coeff( n , k );
+template <typename OutputIterator>
+inline OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t n,
+                                           const binomial_coeff_table& binomial_coeff,
+                                           OutputIterator out) {
+	--n;
 
-    }
+	for (index_t k = dim + 1; k > 0; --k) {
 
-    return out;
-}
+#ifdef USE_BINARY_SEARCH
+		if (binomial_coeff(n, k) > idx) {
+			index_t count;
 
-std::vector<index_t> vertices_of_simplex(const index_t simplex_index, const index_t dim, const index_t n, const binomial_coeff_table& binomial_coeff) {
-    std::vector<index_t> vertices;
-    get_simplex_vertices( simplex_index, dim, n, binomial_coeff, std::back_inserter(vertices) );
-    return vertices;
-}
+#ifdef USE_EXPONENTIAL_SEARCH
+			for (count = 1; (binomial_coeff(n - count, k) > idx); count = std::min(count << 1, n))
+				;
+#else
+			count = n;
+#endif
 
+			while (count > 0) {
+				index_t i = n;
+				index_t step = count >> 1;
+				i -= step;
+				if (binomial_coeff(i, k) > idx) {
+					n = --i;
+					count -= step + 1;
+				} else
+					count = step;
+			}
+		}
+#else
+		while (binomial_coeff(n, k) > idx) --n;
+#endif
 
+		assert(binomial_coeff(n, k) <= idx);
+		assert(binomial_coeff(n + 1, k) > idx);
 
-#ifdef USE_COEFFICIENTS
-struct entry_t {
-    index_t index;
-    coefficient_t coefficient;
-    
-    entry_t(index_t _index, coefficient_t _coefficient) :
-        index(_index), coefficient(_coefficient) {}
-
-    entry_t(index_t _index) :
-        index(_index), coefficient(1) {}
-    
-    entry_t() :
-    index(0), coefficient(1) {}
-};
+		*out++ = n;
 
-inline entry_t make_entry(index_t _index, coefficient_t _coefficient) {
-    return entry_t(_index, _coefficient);
-}
+		idx -= binomial_coeff(n, k);
+	}
 
-inline index_t get_index(entry_t e) {
-    return e.index;
+	return out;
 }
 
-inline index_t get_coefficient(entry_t e) {
-    return e.coefficient;
+std::vector<index_t> vertices_of_simplex(const index_t simplex_index, const index_t dim,
+                                         const index_t n,
+                                         const binomial_coeff_table& binomial_coeff) {
+	std::vector<index_t> vertices;
+	get_simplex_vertices(simplex_index, dim, n, binomial_coeff, std::back_inserter(vertices));
+	return vertices;
 }
 
-inline void set_coefficient(entry_t& e, const coefficient_t c) { e.coefficient = c; }
+#ifdef USE_COEFFICIENTS
+struct entry_t {
+	index_t index;
+	coefficient_t coefficient;
+	entry_t(index_t _index, coefficient_t _coefficient)
+	    : index(_index), coefficient(_coefficient) {}
+	entry_t(index_t _index) : index(_index), coefficient(1) {}
+	entry_t() : index(0), coefficient(1) {}
+};
 
+inline entry_t make_entry(index_t _index, coefficient_t _coefficient) {
+	return entry_t(_index, _coefficient);
+}
+inline index_t get_index(entry_t e) { return e.index; }
+inline index_t get_coefficient(entry_t e) { return e.coefficient; }
+inline void set_coefficient(entry_t& e, const coefficient_t c) { e.coefficient = c; }
 
-bool operator== (const entry_t& e1, const entry_t& e2) {
-    return get_index(e1) == get_index(e2) && get_coefficient(e1) == get_coefficient(e2);
+bool operator==(const entry_t& e1, const entry_t& e2) {
+	return get_index(e1) == get_index(e2) && get_coefficient(e1) == get_coefficient(e2);
 }
 
-std::ostream& operator<< (std::ostream& stream, const entry_t& e) {
-    stream << get_index(e) << ":" << get_coefficient(e);
-    return stream;
+std::ostream& operator<<(std::ostream& stream, const entry_t& e) {
+	stream << get_index(e) << ":" << get_coefficient(e);
+	return stream;
 }
 
 #else
 
 typedef index_t entry_t;
-
-inline const index_t get_index(entry_t i) {
-    return i;
-}
-
-inline index_t get_coefficient(entry_t i) {
-    return 1;
-}
-
-inline entry_t make_entry(index_t _index, coefficient_t _value) {
-    return entry_t(_index);
-}
-
+inline const index_t get_index(entry_t i) { return i; }
+inline index_t get_coefficient(entry_t i) { return 1; }
+inline entry_t make_entry(index_t _index, coefficient_t _value) { return entry_t(_index); }
 inline void set_coefficient(index_t& e, const coefficient_t c) { e = c; }
 
 #endif
 
 inline const entry_t& get_entry(const entry_t& e) { return e; }
 
-template <typename Entry>
-struct greater_index {
-    bool operator() (const Entry& a, const Entry& b) {
-        return get_index(a) > get_index(b);
-    }
+template <typename Entry> struct smaller_index {
+	bool operator()(const Entry& a, const Entry& b) { return get_index(a) < get_index(b); }
 };
 
-
-#ifdef STORE_DIAMETERS
 typedef std::pair<value_t, index_t> diameter_index_t;
-inline value_t get_diameter(diameter_index_t i) {
-    return i.first;
-}
-inline index_t get_index(diameter_index_t i) {
-    return i.second;
-}
-#else
-typedef index_t diameter_index_t;
-#endif
+inline value_t get_diameter(diameter_index_t i) { return i.first; }
+inline index_t get_index(diameter_index_t i) { return i.second; }
 
-#ifdef STORE_DIAMETERS
-class diameter_entry_t: public std::pair<value_t, entry_t> {
+class diameter_entry_t : public std::pair<value_t, entry_t> {
 public:
-    diameter_entry_t(std::pair<value_t, entry_t> p) : std::pair<value_t, entry_t>(p) {}
-#ifdef USE_COEFFICIENTS
-    diameter_entry_t(entry_t e) : std::pair<value_t, entry_t>(0, e) {}
-#endif
-    diameter_entry_t(index_t i) : std::pair<value_t, entry_t>(0, i) {}
+	diameter_entry_t(std::pair<value_t, entry_t> p) : std::pair<value_t, entry_t>(p) {}
+	diameter_entry_t(entry_t e) : std::pair<value_t, entry_t>(0, e) {}
+	diameter_entry_t() : diameter_entry_t(0) {}
 };
 
 inline const entry_t& get_entry(const diameter_entry_t& p) { return p.second; }
 inline entry_t& get_entry(diameter_entry_t& p) { return p.second; }
 inline const index_t get_index(const diameter_entry_t& p) { return get_index(get_entry(p)); }
-inline const coefficient_t get_coefficient(const diameter_entry_t& p) { return get_coefficient(get_entry(p)); }
+inline const coefficient_t get_coefficient(const diameter_entry_t& p) {
+	return get_coefficient(get_entry(p));
+}
 inline const value_t& get_diameter(const diameter_entry_t& p) { return p.first; }
-inline void set_coefficient(diameter_entry_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); }
-inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, coefficient_t _coefficient) {
-    return std::make_pair(_diameter, make_entry(_index, _coefficient));
+inline void set_coefficient(diameter_entry_t& p, const coefficient_t c) {
+	set_coefficient(get_entry(p), c);
 }
-inline diameter_entry_t make_diameter_entry(diameter_index_t _diameter_index, coefficient_t _coefficient) {
-    return std::make_pair(get_diameter(_diameter_index), make_entry(get_index(_diameter_index), _coefficient));
+inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index,
+                                            coefficient_t _coefficient) {
+	return std::make_pair(_diameter, make_entry(_index, _coefficient));
 }
-
-
-struct greater_diameter_or_index {
-    inline bool operator() (const diameter_entry_t& a, const diameter_entry_t& b) {
-        return (get_diameter(a) > get_diameter(b)) ||
-            ((get_diameter(a) == get_diameter(b)) && (get_index(a) > get_index(b)));
-    }
-};
-#else
-typedef entry_t diameter_entry_t;
-inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, coefficient_t _coefficient) {
-    return make_entry(_index, _coefficient);
+inline diameter_entry_t make_diameter_entry(diameter_index_t _diameter_index,
+                                            coefficient_t _coefficient) {
+	return std::make_pair(get_diameter(_diameter_index),
+	                      make_entry(get_index(_diameter_index), _coefficient));
 }
-inline diameter_entry_t make_diameter_entry(index_t _index, coefficient_t _coefficient) {
-    return make_entry(_index, _coefficient);
-}
-#endif
-
 
+template <typename Entry> struct greater_diameter_or_smaller_index {
+	inline bool operator()(const Entry& a, const Entry& b) {
+		return (get_diameter(a) > get_diameter(b)) ||
+		       ((get_diameter(a) == get_diameter(b)) && (get_index(a) < get_index(b)));
+	}
+};
 
-template <typename DistanceMatrix>
-class rips_filtration_comparator {
+template <typename DistanceMatrix> class rips_filtration_comparator {
 public:
-    const DistanceMatrix& dist;
-    
-    const index_t dim;
+	const DistanceMatrix& dist;
+	const index_t dim;
 
 private:
-    mutable std::vector<index_t> vertices;
-    
-    const binomial_coeff_table& binomial_coeff;
-    
-public:
-    rips_filtration_comparator(
-        const DistanceMatrix& _dist,
-        const index_t _dim,
-        const binomial_coeff_table& _binomial_coeff
-    ): dist(_dist), dim(_dim), vertices(_dim + 1), binomial_coeff(_binomial_coeff) {};
-    
-    inline value_t diameter(const index_t index) const {
-        value_t diam = 0;
-        get_simplex_vertices(index, dim, dist.size(), binomial_coeff, vertices.begin() );
-        
-        for (index_t i = 0; i <= dim; ++i)
-            for (index_t j = 0; j < i; ++j) {
-                diam = std::max(diam, dist(vertices[i], vertices[j]));
-            }
-        return diam;
-    }
-
-    inline bool operator()(const index_t a, const index_t b) const
-    {
-        assert(a < binomial_coeff(dist.size(), dim + 1));
-        assert(b < binomial_coeff(dist.size(), dim + 1));
-        
-        value_t a_diam = 0, b_diam = 0;
-    
-        b_diam = diameter(b);
-        
-        get_simplex_vertices(a, dim, dist.size(), binomial_coeff, vertices.begin() );
-        for (index_t i = 0; i <= dim; ++i)
-            for (index_t j = i + 1; j <= dim; ++j) {
-                a_diam = std::max(a_diam, dist(vertices[i], vertices[j]));
-                if (a_diam > b_diam)
-                    return true;
-            }
-       
-        return (a_diam == b_diam) && (a > b);
-    }
-    
-    template <typename Entry>
-    inline bool operator()(const Entry& a, const Entry& b) const
-    {
-        return operator()(get_index(a), get_index(b));
-    }
+	mutable std::vector<index_t> vertices;
+	const binomial_coeff_table& binomial_coeff;
 
+public:
+	rips_filtration_comparator(const DistanceMatrix& _dist, const index_t _dim,
+	                           const binomial_coeff_table& _binomial_coeff)
+	    : dist(_dist), dim(_dim), vertices(_dim + 1), binomial_coeff(_binomial_coeff){};
+
+	inline value_t diameter(const index_t index) const {
+		value_t diam = 0;
+		get_simplex_vertices(index, dim, dist.size(), binomial_coeff, vertices.begin());
+
+		for (index_t i = 0; i <= dim; ++i)
+			for (index_t j = 0; j < i; ++j) {
+				diam = std::max(diam, dist(vertices[i], vertices[j]));
+			}
+		return diam;
+	}
+
+	inline bool operator()(const index_t a, const index_t b) const {
+		assert(a < binomial_coeff(dist.size(), dim + 1));
+		assert(b < binomial_coeff(dist.size(), dim + 1));
+
+		value_t a_diam = 0, b_diam = diameter(b);
+
+		get_simplex_vertices(a, dim, dist.size(), binomial_coeff, vertices.begin());
+		for (index_t i = 0; i <= dim; ++i)
+			for (index_t j = i + 1; j <= dim; ++j) {
+				a_diam = std::max(a_diam, dist(vertices[i], vertices[j]));
+				if (a_diam > b_diam) return true;
+			}
+
+		return (a_diam == b_diam) && (a < b);
+	}
+
+	template <typename Entry> inline bool operator()(const Entry& a, const Entry& b) const {
+		return operator()(get_index(a), get_index(b));
+	}
 };
 
-
 class simplex_coboundary_enumerator {
 private:
-    index_t idx, modified_idx, dim, n, k;
-    
-    const binomial_coeff_table& binomial_coeff;
-    
+	index_t idx, modified_idx, dim, v, k;
+	const binomial_coeff_table& binomial_coeff;
+
 public:
-    inline simplex_coboundary_enumerator (
-        index_t _idx, index_t _dim, index_t _n,
-        const binomial_coeff_table& _binomial_coeff):
-        idx(_idx), modified_idx(_idx), dim(_dim), k(dim + 1), n(_n - 1), binomial_coeff(_binomial_coeff) {}
-    
-    inline bool has_next()
-    {
-        while ( (k != -1 && n != -1) && (binomial_coeff( n , k ) <= idx) ) {
-            idx -= binomial_coeff( n , k );
-            
-            modified_idx -= binomial_coeff( n , k );
-            modified_idx += binomial_coeff( n , k + 1 );
-            
-            --n;
+	inline simplex_coboundary_enumerator(index_t _idx, index_t _dim, index_t _n,
+	                                     const binomial_coeff_table& _binomial_coeff)
+	    : idx(_idx), modified_idx(_idx), dim(_dim), k(dim + 1), v(_n - 1),
+	      binomial_coeff(_binomial_coeff) {}
+
+	inline bool has_next() {
+		while ((v != -1) && (binomial_coeff(v, k) <= idx)) {
+			idx -= binomial_coeff(v, k);
+			modified_idx += binomial_coeff(v, k + 1) - binomial_coeff(v, k);
+			--v;
             --k;
-        }
-        return k != -1 && n != -1;
-    }
-    
-    inline std::pair<entry_t, index_t> next()
-    {
-        while ( binomial_coeff( n , k ) <= idx ) {
-            idx -= binomial_coeff( n , k );
-            
-            modified_idx -= binomial_coeff( n , k );
-            modified_idx += binomial_coeff( n , k + 1 );
-            
-            --n;
-        }
-        auto result = std::make_pair(make_entry(
-            modified_idx + binomial_coeff( n , k + 1 ),
-            k & 1 ? -1 : 1),
-            n);
-        
-        --n;
-        return result;
-    }
+			assert(k != -1);
+		}
+		return v != -1;
+	}
+
+	inline std::pair<entry_t, index_t> next() {
+		auto result =
+		    std::make_pair(make_entry(modified_idx + binomial_coeff(v, k + 1), k & 1 ? -1 : 1), v);
+		--v;
+		return result;
+	}
 };
 
-template <typename DistanceMatrix>
-std::vector<value_t> get_diameters (
-    const DistanceMatrix& dist, const index_t dim,
-    const std::vector<value_t>& previous_diameters,
-    const binomial_coeff_table& binomial_coeff
-    )
-{
-    index_t n = dist.size();
-    
-    std::vector<value_t> diameters(binomial_coeff(n, dim + 1), 0);
-    
-    std::vector<index_t> coboundary;
-
-    for (index_t simplex = 0; simplex < previous_diameters.size(); ++simplex) {
-        coboundary.clear();
-        
-        #ifdef INDICATE_PROGRESS
-        std::cout << "\033[Kpropagating diameter of simplex " << simplex + 1 << "/" << previous_diameters.size() << std::flush << "\r";
-        #endif
-        
-        simplex_coboundary_enumerator coboundaries(simplex, dim - 1, n, binomial_coeff);
-        while (coboundaries.has_next()) {
-            index_t coface = get_index(coboundaries.next().first);
-            diameters[coface] = std::max( diameters[coface], previous_diameters[simplex]);
-        }
-    }
-    
-    #ifdef INDICATE_PROGRESS
-    std::cout << "\033[K";
-    #endif
-    
-    return diameters;
-}
-
-
-class rips_filtration_diameter_comparator {
-private:
-    const std::vector<value_t>& diameters;
-    
-    const index_t dim;
-
-public:
-    std::vector<index_t> vertices;
-    
-    typedef value_t dist_t;
-    
-    const binomial_coeff_table& binomial_coeff;
-    
+class distance_matrix {
 public:
-    rips_filtration_diameter_comparator(
-        const std::vector<value_t>& _diameters,
-        const index_t _dim,
-        const binomial_coeff_table& _binomial_coeff
-        ):
-        diameters(_diameters), dim(_dim),
-        binomial_coeff(_binomial_coeff) {}
-    
-    inline value_t diameter(const index_t a) const {
-        assert(a < diameters.size());
-        return diameters[a];
-    }
-
-    inline bool operator()(const index_t a, const index_t b) const
-    {
-        assert(a < diameters.size());
-        assert(b < diameters.size());
-        
-        dist_t a_diam = diameters[a], b_diam = diameters[b];
-
-        return ((a_diam > b_diam) || ((a_diam == b_diam) && (a > b)));
-    }
-    
-    template <typename Entry>
-    inline bool operator()(const Entry& a, const Entry& b) const
-    {
-        return operator()(get_index(a), get_index(b));
-    }
+	typedef value_t value_type;
+	std::vector<std::vector<value_t>> distances;
+	inline value_t operator()(const index_t a, const index_t b) const { return distances[a][b]; }
+
+	inline size_t size() const { return distances.size(); }
 };
 
+enum compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR };
 
-class distance_matrix  {
+template <compressed_matrix_layout Layout> class compressed_distance_matrix_adapter {
 public:
-    typedef value_t value_type;
-    std::vector<std::vector<value_t> > distances;
-    inline value_t operator()(const index_t a, const index_t b) const {
-        return distances[a][b];
-    }
-    
-    inline size_t size() const {
-        return distances.size();
-    }
-};
+	typedef value_t value_type;
+	std::vector<value_t>& distances;
+	std::vector<value_t*> rows;
 
+	void init_distances() { distances.resize(size() * (size() - 1) / 2); }
 
-class compressed_upper_distance_matrix_adapter  {
-public:
-    typedef value_t value_type;
-    std::vector<value_t>& distances;
-    std::vector<value_t*> rows;
-    
-    index_t n;
-    
-    void init_distances () {
-        distances.resize(n * (n-1) / 2);
-    }
-        
-    void init_rows () {
-        rows.resize(n);
-        value_t* pointer = &distances[0] - 1;
-        for (index_t i = 0; i < n - 1; ++i) {
-            rows[i] = pointer;
-            pointer += n - i - 2;
-        }
-    }
-    
-    compressed_upper_distance_matrix_adapter(std::vector<value_t>& _distances) :
-        distances(_distances)
-    {
-        n = (1 + std::sqrt(1+ 8 * _distances.size())) / 2;
-        
-        assert( distances.size() == n * (n-1) / 2 );
-        
-        init_rows();
-    }
-    
-    inline value_t operator()(const index_t a, const index_t b) const {
-        if (a < b)
-            return rows[a][b];
-        else if (a > b)
-            return rows[b][a];
-        else
-            return 0;
-    }
-    
-    inline size_t size() const {
-        return n;
-    }
-};
+	void init_rows();
 
+	compressed_distance_matrix_adapter(std::vector<value_t>& _distances)
+	    : distances(_distances), rows((1 + std::sqrt(1 + 8 * _distances.size())) / 2) {
+		assert(distances.size() == size() * (size() - 1) / 2);
+		init_rows();
+	}
 
-class compressed_lower_distance_matrix_adapter  {
-public:
-    typedef value_t value_type;
-    std::vector<value_t>& distances;
-    std::vector<value_t*> rows;
-    
-    index_t n;
-    
-    void init_distances () {
-        distances.resize(n * (n-1) / 2);
-    }
-        
-    void init_rows () {
-        rows.resize(n);
-        value_t* pointer = &distances[0];
-        for (index_t i = 1; i < n; ++i) {
-            rows[i] = pointer;
-            pointer += i;
-        }
-    }
-    
-    compressed_lower_distance_matrix_adapter(
-    std::vector<value_t>& _distances, const index_t _n) :
-        distances(_distances), n(_n)
-    {
-        init_distances();
-        init_rows();
-    }
-
-    compressed_lower_distance_matrix_adapter(std::vector<value_t>& _distances) :
-        distances(_distances)
-    {
-        n = (1 + std::sqrt(1+ 8 * distances.size())) / 2;
-        assert( distances.size() == n * (n-1) / 2 );
-        
-        init_rows();
-    }
-    
-    template <typename DistanceMatrix>
-    compressed_lower_distance_matrix_adapter(
-    std::vector<value_t>& _distances, const DistanceMatrix& mat) :
-        distances(_distances), n(mat.size()) {
-        init_distances();
-        init_rows();
-        
-        for (index_t i = 1; i < n; ++i) {
-            for (index_t j = 0; j < i; ++j) {
-                rows[i][j] = mat(i, j);
-            }
-        }
-    }
-    
-    inline value_t operator()(const index_t i, const index_t j) const {
-        if (i > j)
-            return rows[i][j];
-        else if (i < j)
-            return rows[j][i];
-        else
-            return 0;
-    }
-    
-    inline size_t size() const {
-        return n;
-    }
-};
+	template <typename DistanceMatrix>
+	compressed_distance_matrix_adapter(std::vector<value_t>& _distances, const DistanceMatrix& mat)
+	    : distances(_distances), rows(mat.size()) {
+		init_distances();
+		init_rows();
 
+		for (index_t i = 1; i < size(); ++i) {
+			for (index_t j = 0; j < i; ++j) { rows[i][j] = mat(i, j); }
+		}
+	}
 
+	inline value_t operator()(const index_t i, const index_t j) const;
 
+	inline size_t size() const { return rows.size(); }
+};
 
+template <> void compressed_distance_matrix_adapter<LOWER_TRIANGULAR>::init_rows() {
+	value_t* pointer = &distances[0];
+	for (index_t i = 1; i < size(); ++i) {
+		rows[i] = pointer;
+		pointer += i;
+	}
+}
 
-#ifdef USE_COEFFICIENTS
-template <typename Heap>
-inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus)
-{
-    
-    if( column.empty() )
-        return -1;
-    else {
-        auto pivot = column.top();
-        
-        coefficient_t coefficient = 0;
-        do {
-            coefficient = (coefficient + get_coefficient(column.top())) % modulus;
-            column.pop();
-            
-            if( coefficient == 0 ) {
-                if (column.empty()) {
-                    return -1;
-                } else {
-                    pivot = column.top();
-                }
-            }
-        } while ( !column.empty() && get_index(column.top()) == get_index(pivot) );
-        if( get_index(pivot) != -1 ) {
-            set_coefficient(pivot, coefficient);
-        }
-        return pivot;
-    }
+template <> void compressed_distance_matrix_adapter<UPPER_TRIANGULAR>::init_rows() {
+	value_t* pointer = &distances[0] - 1;
+	for (index_t i = 0; i < size() - 1; ++i) {
+		rows[i] = pointer;
+		pointer += size() - i - 2;
+	}
 }
 
-#else
+template <>
+inline value_t compressed_distance_matrix_adapter<UPPER_TRIANGULAR>::
+operator()(const index_t i, const index_t j) const {
+	if (i < j)
+		return rows[i][j];
+	else if (i > j)
+		return rows[j][i];
+	else
+		return 0;
+}
 
-template <typename Heap>
-inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus)
-{
-    
-    if( column.empty() )
-        return -1;
-    else {
-        auto pivot = column.top();
-        column.pop();
-        while( !column.empty() && get_index(column.top()) == get_index(pivot) ) {
-            column.pop();
-            if( column.empty() )
-                return -1;
-            else {
-                pivot = column.top();
-                column.pop();
-            }
-        }
-        return pivot;
-    }
+template <>
+inline value_t compressed_distance_matrix_adapter<LOWER_TRIANGULAR>::
+operator()(const index_t i, const index_t j) const {
+	if (i > j)
+		return rows[i][j];
+	else if (i < j)
+		return rows[j][i];
+	else
+		return 0;
 }
 
-#endif
+typedef compressed_distance_matrix_adapter<LOWER_TRIANGULAR>
+    compressed_lower_distance_matrix_adapter;
+typedef compressed_distance_matrix_adapter<UPPER_TRIANGULAR>
+    compressed_upper_distance_matrix_adapter;
 
-template <typename Heap>
-inline diameter_entry_t get_pivot(Heap& column, coefficient_t modulus)
-{
-    diameter_entry_t result = pop_pivot(column, modulus);
-    if( get_index(result) != -1 ) {
-        column.push(result);
-    }
-    return result;
+
+template <typename Heap> inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) {
+	if (column.empty())
+		return diameter_entry_t(-1);
+	else {
+		auto pivot = column.top();
+        
+#ifdef USE_COEFFICIENTS
+		coefficient_t coefficient = 0;
+		do {
+			coefficient = (coefficient + get_coefficient(column.top())) % modulus;
+			column.pop();
+
+			if (coefficient == 0) {
+				if (column.empty()) {
+					return diameter_entry_t(-1);
+				} else {
+					pivot = column.top();
+				}
+			}
+		} while (!column.empty() && get_index(column.top()) == get_index(pivot));
+		if (get_index(pivot) != -1) { set_coefficient(pivot, coefficient); }
+#else
+		column.pop();
+		while (!column.empty() && get_index(column.top()) == get_index(pivot)) {
+			column.pop();
+			if (column.empty())
+				return diameter_entry_t(-1);
+			else {
+				pivot = column.top();
+				column.pop();
+			}
+		}
+#endif
+		return pivot;
+	}
 }
 
 
+template <typename Heap> inline diameter_entry_t get_pivot(Heap& column, coefficient_t modulus) {
+	diameter_entry_t result = pop_pivot(column, modulus);
+	if (get_index(result) != -1) { column.push(result); }
+	return result;
+}
 
 template <typename Comparator>
-void assemble_columns_to_reduce (
-    std::vector<diameter_index_t>& columns_to_reduce,
-    google::sparse_hash_map<index_t, index_t>& pivot_column_index,
-    const Comparator& comp,
-    index_t dim, index_t n,
-    value_t threshold,
-    const binomial_coeff_table& binomial_coeff
-) {
-    index_t num_simplices = binomial_coeff(n, dim + 2);
-    
-    columns_to_reduce.clear();
-    
-    for (index_t index = 0; index < num_simplices; ++index) {
-        if (pivot_column_index.find(index) == pivot_column_index.end()) {
-            value_t diameter = comp.diameter(index);
-            if (diameter <= threshold)
-                #ifdef STORE_DIAMETERS
-                columns_to_reduce.push_back(std::make_pair(diameter, index));
-                #else
-                columns_to_reduce.push_back(index);
-                #endif
-
-        }
-    }
-
-    #ifdef STORE_DIAMETERS
-    std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), std::greater<diameter_index_t>());
-    #else
-    std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), comp);
-    #endif
+void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce,
+                                hash_map<index_t, index_t>& pivot_column_index,
+                                const Comparator& comp, index_t dim, index_t n, value_t threshold,
+                                const binomial_coeff_table& binomial_coeff) {
+	index_t num_simplices = binomial_coeff(n, dim + 2);
+
+	columns_to_reduce.clear();
+
+	for (index_t index = 0; index < num_simplices; ++index) {
+		if (pivot_column_index.find(index) == pivot_column_index.end()) {
+			value_t diameter = comp.diameter(index);
+			if (diameter <= threshold)
+				columns_to_reduce.push_back(std::make_pair(diameter, index));
+		}
+	}
+
+	std::sort(columns_to_reduce.begin(), columns_to_reduce.end(),
+	          greater_diameter_or_smaller_index<diameter_index_t>());
 }
 
-
-template <typename ValueType>
-class compressed_sparse_matrix  {
+template <typename ValueType> class compressed_sparse_matrix {
 public:
-    std::vector<size_t> bounds;
-    std::vector<ValueType> entries;
-    
-    
-    inline size_t size() const {
-        return bounds.size();
-    }
-    
-    inline typename std::vector<ValueType>::const_iterator cbegin(size_t index) const {
-        assert(index < size());
-        return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1];
-    }
-
-    inline typename std::vector<ValueType>::const_iterator cend(size_t index) const {
-        assert(index < size());
-        return entries.cbegin() + bounds[index];
-    }
-    
-    template <typename Iterator>
-    inline void append(Iterator begin, Iterator end) {
-        for (Iterator it = begin; it != end; ++it) {
-            entries.push_back(*it);
-        }
-        bounds.push_back(entries.size());
-    }
-
-    inline void append() {
-        bounds.push_back(entries.size());
-    }
-
-    inline void push_back(ValueType e) {
-        assert(0 < size());
-        entries.push_back(e);
-        ++bounds.back();
-    }
-    
-    inline void pop_back() {
-        assert(0 < size());
-        entries.pop_back();
-        --bounds.back();
-    }
-    
-    template <typename Collection>
-    inline void append(const Collection collection) {
-        append(collection.cbegin(), collection.cend());
-    }
-
+	std::vector<size_t> bounds;
+	std::vector<ValueType> entries;
+
+	inline size_t size() const { return bounds.size(); }
+
+	inline typename std::vector<ValueType>::const_iterator cbegin(size_t index) const {
+		assert(index < size());
+		return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1];
+	}
+
+	inline typename std::vector<ValueType>::const_iterator cend(size_t index) const {
+		assert(index < size());
+		return entries.cbegin() + bounds[index];
+	}
+
+	template <typename Iterator> inline void append(Iterator begin, Iterator end) {
+		for (Iterator it = begin; it != end; ++it) { entries.push_back(*it); }
+		bounds.push_back(entries.size());
+	}
+
+	inline void append() { bounds.push_back(entries.size()); }
+
+	inline void push_back(ValueType e) {
+		assert(0 < size());
+		entries.push_back(e);
+		++bounds.back();
+	}
+
+	inline void pop_back() {
+		assert(0 < size());
+		entries.pop_back();
+		--bounds.back();
+	}
+
+	template <typename Collection> inline void append(const Collection collection) {
+		append(collection.cbegin(), collection.cend());
+	}
 };
 
-
 template <typename Heap>
-inline std::vector<entry_t> get_column_vector(Heap column, coefficient_t modulus)
-{
-    std::vector<entry_t> temp_col;
-    entry_t pivot = pop_pivot( column, modulus );
-    while( get_index(pivot) != -1 ) {
-        temp_col.push_back( pivot );
-        pivot = pop_pivot( column, modulus );
-    }
-    return temp_col;
+inline void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) {
+	entry_t e = make_entry(i, c);
+	column.push(std::make_pair(diameter, e));
 }
 
+template <typename DistanceMatrix, typename ComparatorCofaces, typename Comparator>
+void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
+                   hash_map<index_t, index_t>& pivot_column_index,
+                   const DistanceMatrix& dist, const ComparatorCofaces& comp,
+                   const Comparator& comp_prev, index_t dim, index_t n, value_t threshold,
+                   coefficient_t modulus, const std::vector<coefficient_t>& multiplicative_inverse,
+                   const binomial_coeff_table& binomial_coeff) {
+
+#ifdef PRINT_PERSISTENCE_PAIRS
+	std::cout << "persistence intervals in dim " << dim << ":" << std::endl;
+#endif
 
-template <typename Heap>
-inline std::vector<entry_t> get_heap_vector(Heap heap)
-{
-    std::vector<entry_t> temp_col;
-    while( !heap.empty() ) {
-        temp_col.push_back( heap.top() );
-        heap.pop();
-    }
-    return temp_col;
-}
+#ifdef ASSEMBLE_REDUCTION_MATRIX
+	compressed_sparse_matrix<diameter_entry_t> reduction_matrix;
+#else
+#ifdef USE_COEFFICIENTS
+	std::vector<diameter_entry_t> reduction_coefficients;
+#endif
+#endif
 
+    std::vector<diameter_entry_t> coface_entries;
+    std::vector<index_t> vertices;
 
-template <typename Heap>
-inline void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) {
-    entry_t e = make_entry(i, c);
-#ifdef STORE_DIAMETERS
-    column.push(std::make_pair(diameter, e));
-#else
-    column.push(e);
+	for (index_t i = 0; i < columns_to_reduce.size(); ++i) {
+		auto column_to_reduce = columns_to_reduce[i];
+
+#ifdef ASSEMBLE_REDUCTION_MATRIX
+		std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>,
+		                    smaller_index<diameter_entry_t>> reduction_column;
 #endif
-}
 
+		std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>,
+		                    greater_diameter_or_smaller_index<diameter_entry_t>> working_coboundary;
 
-template <typename DistanceMatrix, typename ComparatorCofaces, typename Comparator>
-void compute_pairs(
-    std::vector<diameter_index_t>& columns_to_reduce,
-    google::sparse_hash_map<index_t, index_t>& pivot_column_index,
-    const DistanceMatrix& dist,
-    const ComparatorCofaces& comp, const Comparator& comp_prev,
-    index_t dim, index_t n,
-    value_t threshold,
-    coefficient_t modulus, const std::vector<coefficient_t>& multiplicative_inverse,
-    const binomial_coeff_table& binomial_coeff
-) {
-    
-    #ifdef PRINT_PERSISTENCE_PAIRS
-    std::cout << "persistence intervals in dim " << dim << ":" << std::endl;
-    #endif
-    
-    #ifdef ASSEMBLE_REDUCTION_MATRIX
-    compressed_sparse_matrix <diameter_entry_t> reduction_matrix;
-    #else
-    #ifdef USE_COEFFICIENTS
-    std::vector <diameter_entry_t> reduction_coefficients;
-    #endif
-    #endif
-    
-//    size_t boundary_additions = 0;
-
-    for (index_t i = 0; i < columns_to_reduce.size(); ++i) {
-        auto column_to_reduce = columns_to_reduce[i];
-        
-        #ifdef ASSEMBLE_REDUCTION_MATRIX
-        
-        std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, greater_index<diameter_entry_t>> reduction_column;
-        
-        #endif
-        
-        
-        #ifdef STORE_DIAMETERS
-        std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, greater_diameter_or_index > working_coboundary;
-        #else
-        std::priority_queue<entry_t, std::vector<entry_t>, decltype(comp) > working_coboundary(comp);
-        #endif
-        
-        #ifdef STORE_DIAMETERS
-        value_t diameter = get_diameter(columns_to_reduce[i]);
-        #else
-        value_t diameter = comp_prev.diameter(get_index(column_to_reduce));
-        #endif
+		value_t diameter = get_diameter(columns_to_reduce[i]);
 
-        
-        #ifdef INDICATE_PROGRESS
-        std::cout << "\033[K" << "reducing column " << i + 1 << "/" << columns_to_reduce.size()
-        << " (diameter " << diameter << ")"
-        << std::flush << "\r";
-        #endif
-        
-        index_t j = i;
-        
-        auto column_to_add = column_to_reduce;
-        
-        
-        // start with a dummy pivot entry with coefficient -1 in order to initialize working_coboundary
-        // with the coboundary of the simplex with index column_to_reduce
-        
-        diameter_entry_t pivot = make_diameter_entry(0, -1, modulus - 1);
-        
-//        std::cout << "reducing " << column_to_reduce << ":" << std::endl;
-        
-        #ifdef ASSEMBLE_REDUCTION_MATRIX
-        reduction_matrix.append();
-        #endif
-        
+#ifdef INDICATE_PROGRESS
+		std::cout << "\033[K"
+		          << "reducing column " << i + 1 << "/" << columns_to_reduce.size() << " (diameter "
+		          << diameter << ")" << std::flush << "\r";
+#endif
 
-        // initialize reduction_matrix as identity matrix
-        #ifdef ASSEMBLE_REDUCTION_MATRIX
-        reduction_matrix.push_back(make_diameter_entry(column_to_reduce, 1));
-        #else
-        #ifdef USE_COEFFICIENTS
-        reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, 1));
-        #endif
-        #endif
-        
-//        --boundary_additions;
-        
-        do {
-            
-            const coefficient_t factor = modulus - get_coefficient(pivot);
+		index_t j = i;
+		auto column_to_add = column_to_reduce;
 
+		// start with a dummy pivot entry with coefficient -1 in order to initialize
+		// working_coboundary with the coboundary of the simplex with index column_to_reduce
+		diameter_entry_t pivot = make_diameter_entry(0, -1, modulus - 1);
 
-//            std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, decltype(comp) > eliminating_coboundary(comp);
+#ifdef ASSEMBLE_REDUCTION_MATRIX
+		reduction_matrix.append();
+#endif
 
-//            std::cout << "w:" << get_column_vector(working_coboundary, modulus) << std::endl;
+#ifdef ASSEMBLE_REDUCTION_MATRIX
+		// initialize reduction_matrix as identity matrix
+		reduction_matrix.push_back(make_diameter_entry(column_to_reduce, 1));
+#else
+#ifdef USE_COEFFICIENTS
+		reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, 1));
+#endif
+#endif
 
-            #ifdef ASSEMBLE_REDUCTION_MATRIX
-            for (auto it = reduction_matrix.cbegin(j); it != reduction_matrix.cend(j); ++it)
-            #endif
-            {
-                
-//                ++boundary_additions;
-
-
-                #ifdef ASSEMBLE_REDUCTION_MATRIX
-                auto simplex = *it;
-                coefficient_t simplex_coefficient = get_coefficient(simplex);
-                simplex_coefficient *= factor;
-                simplex_coefficient %= modulus;
-                #else
-                #ifdef USE_COEFFICIENTS
-                const auto& simplex = reduction_coefficients[j];
-                #else
-                const auto simplex = column_to_add;
-                #endif
-                #endif
-                
-                #ifdef STORE_DIAMETERS
-                value_t simplex_diameter = get_diameter(simplex);
-                #else
-                value_t simplex_diameter = comp_prev.diameter(get_index(simplex));
-                #endif
-                assert(simplex_diameter == comp_prev.diameter(get_index(simplex)));
-
-                #ifdef ASSEMBLE_REDUCTION_MATRIX
-                reduction_column.push( make_diameter_entry( simplex_diameter, get_index(simplex), simplex_coefficient ) );
-                #endif
-                
-                std::vector<index_t> vertices = vertices_of_simplex(get_index(simplex), dim, n, binomial_coeff);
+		bool might_be_apparent_pair = true;
+        
+		do {
+			const coefficient_t factor = modulus - get_coefficient(pivot);
+
+#ifdef ASSEMBLE_REDUCTION_MATRIX
+			for (auto it = reduction_matrix.cbegin(j); it != reduction_matrix.cend(j); ++it)
+#endif
+			{
+#ifdef ASSEMBLE_REDUCTION_MATRIX
+				auto simplex = *it;
+				coefficient_t simplex_coefficient = get_coefficient(simplex);
+				simplex_coefficient *= factor;
+				simplex_coefficient %= modulus;
+#else
+#ifdef USE_COEFFICIENTS
+				const auto& simplex = reduction_coefficients[j];
+#else
+				const auto simplex = column_to_add;
+#endif
+#endif
+				value_t simplex_diameter = get_diameter(simplex);
+				assert(simplex_diameter == comp_prev.diameter(get_index(simplex)));
+
+#ifdef ASSEMBLE_REDUCTION_MATRIX
+				reduction_column.push(
+				    make_diameter_entry(simplex_diameter, get_index(simplex), simplex_coefficient));
+#endif
+
+                vertices.clear();
+				get_simplex_vertices(get_index(simplex), dim, n, binomial_coeff, std::back_inserter(vertices));
                 
-                simplex_coboundary_enumerator cofaces(get_index(simplex), dim, n, binomial_coeff);
-                while (cofaces.has_next()) {
-                    auto coface_descriptor = cofaces.next();
-                    
-                    entry_t coface = coface_descriptor.first;
-                    index_t covertex = coface_descriptor.second;
-
-                    index_t coface_index = get_index(coface);
-
-                    value_t coface_diameter = simplex_diameter;
-                    for (index_t v : vertices) {
-                        coface_diameter = std::max(coface_diameter, dist(v, covertex));
-                    }
-                    
-                    assert(comp.diameter(coface_index) == coface_diameter);
-                    
-                    if (coface_diameter <= threshold) {
-                        coefficient_t coface_coefficient = get_coefficient(coface) * get_coefficient(simplex) * factor;
-                        coface_coefficient %= modulus;
-                        if (coface_coefficient < 0) coface_coefficient += modulus;
+                coface_entries.clear();
+				simplex_coboundary_enumerator cofaces(get_index(simplex), dim, n, binomial_coeff);
+				while (cofaces.has_next()) {
+					auto coface_descriptor = cofaces.next();
+					entry_t coface = coface_descriptor.first;
+					index_t covertex = coface_descriptor.second;
+					index_t coface_index = get_index(coface);
+					value_t coface_diameter = simplex_diameter;
+					for (index_t v : vertices) {
+						coface_diameter = std::max(coface_diameter, dist(v, covertex));
+					}
+					assert(comp.diameter(coface_index) == coface_diameter);
+
+					if (coface_diameter <= threshold) {
+						coefficient_t coface_coefficient =
+						    get_coefficient(coface) * get_coefficient(simplex) * factor;
+						coface_coefficient %= modulus;
+						if (coface_coefficient < 0) coface_coefficient += modulus;
+						assert(coface_coefficient >= 0);
                         
-                        assert(coface_coefficient >= 0);
+                        diameter_entry_t coface_entry = make_diameter_entry(coface_diameter, coface_index, coface_coefficient);
+                        coface_entries.push_back(coface_entry);
                         
-                        push_entry(working_coboundary, coface_index, coface_coefficient, coface_diameter);
-//                        push_entry(eliminating_coboundary, coface_index, coface_coefficient, coface_diameter);
-                    }
-                }
-            }
+                        if (might_be_apparent_pair && (simplex_diameter == coface_diameter)) {
+                            if (pivot_column_index.find(coface_index) == pivot_column_index.end())
+                            {
+                                pivot = coface_entry;
+                                goto found_pivot;
+                            }
+                            might_be_apparent_pair = false;
+                        }
+					}
+				}
+                for (auto e: coface_entries) working_coboundary.push(e);
+			}
             
+			pivot = get_pivot(working_coboundary, modulus);
 
+        found_pivot:
+			if (get_index(pivot) != -1) {
+				auto pair = pivot_column_index.find(get_index(pivot));
 
-            
-//            std::cout << get_heap_vector(working_coboundary) << std::endl;
+				if (pair == pivot_column_index.end()) {
+					pivot_column_index.insert(std::make_pair(get_index(pivot), i));
 
-//            std::cout << "e:" << get_column_vector(eliminating_coboundary, modulus) << std::endl;
-//            std::cout << "w:" << get_column_vector(working_coboundary, modulus) << std::endl << std::endl;
-            
-            pivot = get_pivot(working_coboundary, modulus);
-            
-//            std::cout << "pivot " << get_index(pivot) << std::endl;
-        
-            if (get_index(pivot) != -1) {
-                auto pair = pivot_column_index.find(get_index(pivot));
-                
-                if (pair == pivot_column_index.end()) {
-//                    std::cout << std::endl;
-                    
-                    pivot_column_index.insert(std::make_pair(get_index(pivot), i));
-                    
-                    #ifdef USE_COEFFICIENTS
-                    const coefficient_t inverse = multiplicative_inverse[ get_coefficient( pivot ) ];
-                    #else
-                    #ifdef ASSEMBLE_REDUCTION_MATRIX
-                    const coefficient_t inverse = 1;
-                    #endif
-                    #endif
-                    
-                    // replace current column of reduction_matrix (with a single diagonal 1 entry)
-                    // by reduction_column (possibly with a different entry on the diagonal)
-                    #ifdef ASSEMBLE_REDUCTION_MATRIX
-                    reduction_matrix.pop_back();
-                    while (true) {
-                        diameter_entry_t e = pop_pivot(reduction_column, modulus);
-                        index_t index = get_index(e);
-                        if (index == -1)
-                            break;
-                        const coefficient_t coefficient = inverse * get_coefficient(e) % modulus;
-                        assert(coefficient > 0);
-                        reduction_matrix.push_back(make_diameter_entry(get_diameter(e), index, coefficient));
-                    }
-                    #else
-                    #ifdef USE_COEFFICIENTS
-                    reduction_coefficients.pop_back();
-                    reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, inverse));
-                    #endif
-                    #endif
-                    
-                    #ifdef PRINT_PERSISTENCE_PAIRS
-                    #ifdef STORE_DIAMETERS
-                    value_t death = get_diameter(pivot);
-                    #else
-                    value_t death = comp.diameter(get_index(pivot));
-                    #endif
-                    if (diameter != death)
-                    {
-                        #ifdef INDICATE_PROGRESS
-                        std::cout << "\033[K";
-                        #endif
-                        std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush;
-                        
-//                        std::cout << " (" << vertices_of_simplex(column_to_reduce, dim, n, binomial_coeff) << ", " << vertices_of_simplex(get_index(pivot), dim + 1, n, binomial_coeff) << ")" << std::endl;
-                    }
-                    #endif
+#ifdef USE_COEFFICIENTS
+					const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)];
+#endif
 
-                    break;
-                }
-                
-                j = pair->second;
-                column_to_add = columns_to_reduce[j];
-            }
+#ifdef ASSEMBLE_REDUCTION_MATRIX
+					// replace current column of reduction_matrix (with a single diagonal 1 entry)
+					// by reduction_column (possibly with a different entry on the diagonal)
+					reduction_matrix.pop_back();
+					while (true) {
+						diameter_entry_t e = pop_pivot(reduction_column, modulus);
+						index_t index = get_index(e);
+						if (index == -1) break;
+#ifdef USE_COEFFICIENTS
+						const coefficient_t coefficient = inverse * get_coefficient(e) % modulus;
+						assert(coefficient > 0);
+#else
+						const coefficient_t coefficient = 1;
+#endif
+						reduction_matrix.push_back(
+						    make_diameter_entry(get_diameter(e), index, coefficient));
+					}
+#else
+#ifdef USE_COEFFICIENTS
+					reduction_coefficients.pop_back();
+					reduction_coefficients.push_back(
+					    make_diameter_entry(column_to_reduce, inverse));
+#endif
+#endif
 
-        } while ( get_index(pivot) != -1 );
-        
-        #ifdef PRINT_PERSISTENCE_PAIRS
-        if ( get_index(pivot) == -1 ) {
-//            std::cout << std::endl;
-        
-            value_t birth = diameter;
-            #ifdef INDICATE_PROGRESS
-            std::cout << "\033[K";
-            #endif
-            std::cout << " [" << birth << ", )" << std::endl << std::flush;
-        }
-        #endif
-
-    
-//        #ifdef ASSEMBLE_REDUCTION_MATRIX
-//        std::cout << "reduction matrix fill-in: " << i + 1 << " + " << reduction_matrix.entries.size() - (i + 1) << std::endl;
-//        #endif
-    
-    }
-
-    #ifdef INDICATE_PROGRESS
-    std::cout << "\033[K";
-    #endif
-    
-//    std::cout << boundary_additions << " boundary additions" << std::endl;
+#ifdef PRINT_PERSISTENCE_PAIRS
+					value_t death = get_diameter(pivot);
+					if (diameter != death) {
+#ifdef INDICATE_PROGRESS
+						std::cout << "\033[K";
+#endif
+						std::cout << " [" << diameter << "," << death << ")" << std::endl
+						          << std::flush;
+					}
+#endif
+					break;
+				}
+				j = pair->second;
+				column_to_add = columns_to_reduce[j];
+			}
+		} while (get_index(pivot) != -1);
+
+#ifdef PRINT_PERSISTENCE_PAIRS
+		if (get_index(pivot) == -1) {
+			value_t birth = diameter;
+#ifdef INDICATE_PROGRESS
+			std::cout << "\033[K";
+#endif
+			std::cout << " [" << birth << ", )" << std::endl << std::flush;
+		}
+#endif
+	}
 
+#ifdef INDICATE_PROGRESS
+	std::cout << "\033[K";
+#endif
 }
 
-bool is_prime(const long n) {
-    bool is_prime = true;
-    for (int i = 2; i <= n/2; ++i)
-        if (n%i == 0) {
-            is_prime = false;
-            break;
-        }
-    return is_prime;
+bool is_prime(const coefficient_t n) {
+	if (!(n & 1) || n < 2) return n == 2;
+	for (coefficient_t p = 3, q = n / p, r = n % p; p <= q; p += 2, q = n / p, r = n % p)
+		if (!r) return false;
+	return true;
 }
 
-void print_usage_and_exit(int exit_code)
-{
-    std::cerr << "Usage: " << "ripser " << "[options] filename" << std::endl;
-    std::cerr << std::endl;
-    std::cerr << "Options:" << std::endl;
-    std::cerr << std::endl;
-    std::cerr << "  --help           print this screen" << std::endl;
-    std::cerr << "  --top_dim <k>    compute persistent homology up to dimension <k>" << std::endl;
-    std::cerr << "  --threshold <t>  compute Rips complexes up to diameter <t>" << std::endl;
-    #ifdef USE_COEFFICIENTS
-    std::cerr << "  --modulus <p>    compute homology with coefficients in the prime field Z/<p>Z" << std::endl;
-    #endif
-    
-    exit(exit_code);
-}
+void print_usage_and_exit(int exit_code) {
+	std::cerr << "Usage: "
+	          << "ripser "
+	          << "[options] filename" << std::endl;
+	std::cerr << std::endl;
+	std::cerr << "Options:" << std::endl;
+	std::cerr << std::endl;
+	std::cerr << "  --help           print this screen" << std::endl;
+	std::cerr << "  --top_dim <k>    compute persistent homology up to dimension <k>" << std::endl;
+	std::cerr << "  --threshold <t>  compute Rips complexes up to diameter <t>" << std::endl;
+#ifdef USE_COEFFICIENTS
+	std::cerr << "  --modulus <p>    compute homology with coefficients in the prime field Z/<p>Z"
+	          << std::endl;
+#endif
 
-int main( int argc, char** argv ) {
-    
-    if( argc < 2 ) print_usage_and_exit(-1);
-
-    const char *filename = nullptr;
-    
-    index_t dim_max = 1;
-    value_t threshold = std::numeric_limits<value_t>::max();
-
-    #ifdef USE_COEFFICIENTS
-    coefficient_t modulus = 2;
-    #else
-    const coefficient_t modulus = 2;
-    #endif
-
-    for( index_t i = 1; i < argc; ++i ) {
-        const std::string arg(argv[ i ]);
-        if( arg == "--help" ) {
-            print_usage_and_exit(0);
-        } else if( arg == "--top_dim" ) {
-            std::string parameter = std::string( argv[ ++i ] );
-            size_t next_pos;
-            dim_max = std::stol( parameter, &next_pos );
-            if( next_pos != parameter.size() )
-                print_usage_and_exit( -1 );
-        } else if( arg == "--threshold" ) {
-            std::string parameter = std::string( argv[ ++i ] );
-            size_t next_pos;
-            threshold = std::stof( parameter, &next_pos );
-            if( next_pos != parameter.size() )
-                print_usage_and_exit( -1 );
-        #ifdef USE_COEFFICIENTS
-        } else if( arg == "--modulus" ) {
-            std::string parameter = std::string( argv[ ++i ] );
-            size_t next_pos;
-            modulus = std::stol( parameter, &next_pos );
-            if( next_pos != parameter.size() || !is_prime(modulus) )
-                print_usage_and_exit( -1 );
-        #endif
-        } else {
-            if (filename) {
-                print_usage_and_exit( -1 );
-            }
-            filename = argv[i];
-        }
-    }
-
-
-    std::ifstream input_stream( filename, std::ios_base::binary | std::ios_base::in );
-    if( input_stream.fail( ) ) {
-        std::cerr << "couldn't open file " << filename << std::endl;
-        exit(-1);
-    }
-    
-    std::vector<std::vector<value_t>> diameters(dim_max + 2);
-    
-    #ifdef FILE_FORMAT_DIPHA
-    
-    int64_t magic_number;
-    input_stream.read( reinterpret_cast<char*>(&magic_number), sizeof( int64_t ) );
-    if( magic_number != 8067171840 ) {
-        std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl;
-        exit(-1);
-    }
-
-    int64_t file_type;
-    input_stream.read( reinterpret_cast<char*>(&file_type), sizeof( int64_t ) );
-    if( file_type != 7 ) {
-        std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl;
-        exit(-1);
-    }
-   
-    int64_t n;
-    input_stream.read( reinterpret_cast<char*>(&n), sizeof( int64_t ) );
-    
-    distance_matrix dist_full;
-    dist_full.distances = std::vector<std::vector<value_t>>(n, std::vector<value_t>(n));
-    
-    for( int i = 0; i < n; ++i ) {
-        for( int j = 0; j < n; ++j ) {
-            double val;
-            input_stream.read( reinterpret_cast<char*>(&val), sizeof(double) );
-            dist_full.distances[i][j] = val;
-        }
-    }
-    
-    std::cout << "distance matrix with " << n << " points" << std::endl;
-
-    compressed_lower_distance_matrix_adapter dist(diameters[1], dist_full);
-
-    std::cout << "distance matrix transformed to lower triangular form" << std::endl;
-
-    #endif
-
-
-    #ifdef FILE_FORMAT_UPPER_TRIANGULAR_CSV
-    
-    std::vector<value_t> distances;
-    std::string value_string;
-    while(std::getline(input_stream, value_string, ','))
-        distances.push_back(std::stod(value_string));
-    
-    compressed_upper_distance_matrix_adapter dist_upper(distances);
-    
-    index_t n = dist_upper.size();
-    
-    std::cout << "distance matrix with " << n << " points" << std::endl;
- 
-    compressed_lower_distance_matrix_adapter dist(diameters[1], dist_upper);
-
-    std::cout << "distance matrix transformed to lower triangular form" << std::endl;
-    
-    #endif
-    
-
-    #ifdef FILE_FORMAT_LOWER_TRIANGULAR_CSV
-    
-    std::vector<value_t>& distances = diameters[1];
-    std::string value_string;
-    while(std::getline(input_stream, value_string, ','))
-        distances.push_back(std::stod(value_string));
-    
-    compressed_lower_distance_matrix_adapter dist(distances);
-    
-    index_t n = dist.size();
-    
-    std::cout << "distance matrix with " << n << " points" << std::endl;
-     
-    #endif
-
-
-
-    auto result = std::minmax_element(dist.distances.begin(), dist.distances.end());
-    std::cout << "value range: [" << *result.first << "," << *result.second << "]" << std::endl;
-    
-
-    assert(dim_max <= n - 2);
-    
-    binomial_coeff_table binomial_coeff(n, dim_max + 2);
-    
-    std::vector<coefficient_t> multiplicative_inverse(multiplicative_inverse_vector(modulus));
-    
-    std::vector<diameter_index_t> columns_to_reduce;
-    
-    
-    for (index_t index = n; index-- > 0; ) {
-        #ifdef STORE_DIAMETERS
-        columns_to_reduce.push_back(std::make_pair(0, index));
-        #else
-        columns_to_reduce.push_back(index);
-        #endif
-
-    }
-
-    
-
-    index_t dim = 0;
-
-    {
-        rips_filtration_diameter_comparator comp(diameters[1], dim + 1, binomial_coeff);
-        rips_filtration_comparator<decltype(dist)> comp_prev(dist, dim, binomial_coeff);
-
-        google::sparse_hash_map<index_t, index_t> pivot_column_index;
-
-        compute_pairs(
-            columns_to_reduce,
-            pivot_column_index,
-            dist,
-            comp, comp_prev,
-            dim, n,
-            threshold,
-            modulus, multiplicative_inverse,
-            binomial_coeff
-        );
-        
-        assemble_columns_to_reduce(
-            columns_to_reduce,
-            pivot_column_index,
-            comp,
-            dim, n,
-            threshold,
-            binomial_coeff
-        );
-    }
-    
-    #ifdef PRECOMPUTE_DIAMETERS
-    #ifdef PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION
-    for (dim = 1; dim <= dim_max; ++dim) {
-    #else
-    for (dim = 1; dim < dim_max; ++dim) {
-    #endif
-    #else
-    for (dim = 1; dim <= dim_max; ++dim) {
-    #endif
-    
-        #ifdef PRECOMPUTE_DIAMETERS
-        
-        #ifdef INDICATE_PROGRESS
-        std::cout << "precomputing " << dim + 1 << "-simplex diameters" << std::endl;
-        #endif
-        diameters[dim + 1] = get_diameters( dist, dim + 1, diameters[dim], binomial_coeff );
-    
-        rips_filtration_diameter_comparator comp(diameters[dim + 1], dim + 1, binomial_coeff);
-        rips_filtration_diameter_comparator comp_prev(diameters[dim], dim, binomial_coeff);
-        
-        #else
-        
-        rips_filtration_comparator<decltype(dist)> comp(dist, dim + 1, binomial_coeff);
-        rips_filtration_comparator<decltype(dist)> comp_prev(dist, dim, binomial_coeff);
-        
-        #endif
-
-        google::sparse_hash_map<index_t, index_t> pivot_column_index;
-        pivot_column_index.resize(columns_to_reduce.size());
-    
-        compute_pairs(
-            columns_to_reduce,
-            pivot_column_index,
-            dist,
-            comp, comp_prev,
-            dim, n,
-            threshold,
-            modulus,
-            multiplicative_inverse,
-            binomial_coeff
-        );
-        
-        if (dim < dim_max)
-            assemble_columns_to_reduce(
-                columns_to_reduce,
-                pivot_column_index,
-                comp,
-                dim, n,
-                threshold,
-                binomial_coeff
-            );
-        
-//        if ( dim > 1 )
-//            diameters[dim] = std::vector<value_t>();
-    }
-        
-    #ifdef PRECOMPUTE_DIAMETERS
-    #ifndef PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION
-    {
-        dim = dim_max;
-        
-        rips_filtration_diameter_comparator comp_prev(diameters[dim], dim, binomial_coeff);
-        rips_filtration_comparator<decltype(dist)> comp(dist, dim + 1, binomial_coeff);
-        
-        std::unordered_map<index_t, index_t> pivot_column_index;
-        pivot_column_index.resize(columns_to_reduce.size());
-
-        compute_pairs(
-            columns_to_reduce,
-            pivot_column_index,
-            comp, comp_prev,
-            dim, n,
-            threshold,
-            modulus,
-            multiplicative_inverse,
-            binomial_coeff
-        );
-    }
-    #endif
-    #endif
-    
+	exit(exit_code);
 }
+
+int main(int argc, char** argv) {
+
+	if (argc < 2) print_usage_and_exit(-1);
+
+	const char* filename = nullptr;
+
+	index_t dim_max = 1;
+	value_t threshold = std::numeric_limits<value_t>::max();
+
+#ifdef USE_COEFFICIENTS
+	coefficient_t modulus = 2;
+#else
+	const coefficient_t modulus = 2;
+#endif
+
+	for (index_t i = 1; i < argc; ++i) {
+		const std::string arg(argv[i]);
+		if (arg == "--help") {
+			print_usage_and_exit(0);
+		} else if (arg == "--top_dim") {
+			std::string parameter = std::string(argv[++i]);
+			size_t next_pos;
+			dim_max = std::stol(parameter, &next_pos);
+			if (next_pos != parameter.size()) print_usage_and_exit(-1);
+		} else if (arg == "--threshold") {
+			std::string parameter = std::string(argv[++i]);
+			size_t next_pos;
+			threshold = std::stof(parameter, &next_pos);
+			if (next_pos != parameter.size()) print_usage_and_exit(-1);
+#ifdef USE_COEFFICIENTS
+		} else if (arg == "--modulus") {
+			std::string parameter = std::string(argv[++i]);
+			size_t next_pos;
+			modulus = std::stol(parameter, &next_pos);
+			if (next_pos != parameter.size() || !is_prime(modulus)) print_usage_and_exit(-1);
+#endif
+		} else {
+			if (filename) { print_usage_and_exit(-1); }
+			filename = argv[i];
+		}
+	}
+
+	std::ifstream input_stream(filename, std::ios_base::binary | std::ios_base::in);
+	if (input_stream.fail()) {
+		std::cerr << "couldn't open file " << filename << std::endl;
+		exit(-1);
+	}
+
+	std::vector<value_t> diameters;
+
+#ifdef FILE_FORMAT_DIPHA
+
+	int64_t magic_number;
+	input_stream.read(reinterpret_cast<char*>(&magic_number), sizeof(int64_t));
+	if (magic_number != 8067171840) {
+		std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl;
+		exit(-1);
+	}
+
+	int64_t file_type;
+	input_stream.read(reinterpret_cast<char*>(&file_type), sizeof(int64_t));
+	if (file_type != 7) {
+		std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl;
+		exit(-1);
+	}
+
+	int64_t n;
+	input_stream.read(reinterpret_cast<char*>(&n), sizeof(int64_t));
+
+	distance_matrix dist_full;
+	dist_full.distances = std::vector<std::vector<value_t>>(n, std::vector<value_t>(n));
+
+	for (int i = 0; i < n; ++i) {
+		for (int j = 0; j < n; ++j) {
+			double val;
+			input_stream.read(reinterpret_cast<char*>(&val), sizeof(double));
+			dist_full.distances[i][j] = val;
+		}
+	}
+	std::cout << "distance matrix with " << n << " points" << std::endl;
+
+	compressed_lower_distance_matrix_adapter dist(diameters, dist_full);
+	std::cout << "distance matrix transformed to lower triangular form" << std::endl;
+
+#endif
+
+#ifdef FILE_FORMAT_UPPER_TRIANGULAR_CSV
+
+	std::vector<value_t> distances;
+	std::string value_string;
+	while (std::getline(input_stream, value_string, ','))
+		distances.push_back(std::stod(value_string));
+
+	compressed_upper_distance_matrix_adapter dist_upper(distances);
+
+	index_t n = dist_upper.size();
+	std::cout << "distance matrix with " << n << " points" << std::endl;
+
+	compressed_lower_distance_matrix_adapter dist(diameters, dist_upper);
+	std::cout << "distance matrix transformed to lower triangular form" << std::endl;
+
+#endif
+
+#ifdef FILE_FORMAT_LOWER_TRIANGULAR_CSV
+
+	std::vector<value_t>& distances = diameters;
+	std::string value_string;
+	while (std::getline(input_stream, value_string, ','))
+		distances.push_back(std::stod(value_string));
+
+	compressed_lower_distance_matrix_adapter dist(distances);
+
+	index_t n = dist.size();
+
+	std::cout << "distance matrix with " << n << " points" << std::endl;
+
+#endif
+
+	auto result = std::minmax_element(dist.distances.begin(), dist.distances.end());
+	std::cout << "value range: [" << *result.first << "," << *result.second << "]" << std::endl;
+
+	assert(dim_max <= n - 2);
+
+	binomial_coeff_table binomial_coeff(n, dim_max + 2);
+	std::vector<coefficient_t> multiplicative_inverse(multiplicative_inverse_vector(modulus));
+
+	std::vector<diameter_index_t> columns_to_reduce;
+	for (index_t index = n; index-- > 0;) columns_to_reduce.push_back(diameter_index_t(0, index));
+
+	for (index_t dim = 0; dim <= dim_max; ++dim) {
+
+		rips_filtration_comparator<decltype(dist)> comp(dist, dim + 1, binomial_coeff);
+		rips_filtration_comparator<decltype(dist)> comp_prev(dist, dim, binomial_coeff);
+
+		hash_map<index_t, index_t> pivot_column_index;
+		pivot_column_index.reserve(columns_to_reduce.size());
+
+		compute_pairs(columns_to_reduce, pivot_column_index, dist, comp, comp_prev, dim, n,
+		              threshold, modulus, multiplicative_inverse, binomial_coeff);
+
+		if (dim < dim_max)
+			assemble_columns_to_reduce(columns_to_reduce, pivot_column_index, comp, dim, n,
+			                           threshold, binomial_coeff);
+	}
+}
\ No newline at end of file