+
+/** Counts the number of common elements in two sorted vectors. Equal to
+ * counting the number of elements given by std::set_intersection. */
+template <class InputIterator1, class InputIterator2>
+size_t count_set_intersection (InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, InputIterator2 last2)
+{
+ size_t count = 0;
+
+ // As long as we did not either end, look for common elements
+ while (first1!=last1 && first2!=last2)
+ {
+ if (*first1 < *first2)
+ ++first1;
+ else if (*first2 < *first1)
+ ++first2;
+ // Element is common to both
+ else {
+ ++count;
+ ++first1;
+ ++first2;
+ }
+ }
+ return count;
+}
+
+/** Multiply a*b and save the product in 'result'. It is assumed that 'result' is intially empty and has appropriate size. */
+template<class IT, class D1, class D2, class RT>
+void multiply_matrix(RT &result, const boolean_rowmatrix_base<IT, D1> &a, const boolean_colmatrix_base<IT, D2> &b) {
+ assert(a.height() == b.width());
+
+ for (unsigned r=0; r < a.height(); ++r) {
+ const typename boolean_rowmatrix_base<IT, D1>::row_type &row = a.get_row(r);
+ for (unsigned c=0; c < b.width(); ++c) {
+ const typename boolean_colmatrix_base<IT, D2>::column_type &col = b.get_column(c);
+ if (count_set_intersection(row.begin(), row.end(), col.begin(), col.end()) % 2 == 1)
+ result.set(r, c, true);
+ }
+ }
+}
+
+template<class MT>
+MT create_unit_matrix(size_t size) {
+ MT mat(size);
+ for (unsigned i=0; i < size; ++i)
+ mat.set(i, i, true);
+ return mat;
+}
+