Fix bugs in USE_TWO_CUTS

group-connectivity.h

@@ -64,19 +64,10 @@ struct Tester : public AbstractTester {
     numForb = 1;
     while (edges_-- > 0) numForb *= (Ring::size - 1);
 
-    size_t num_classes = 1;
-    edges_ = spanningTree.size();
-    while (edges_-- > 0) num_classes *= Ring::size;
-    ClassesType C(num_classes, 0);
-    C.swap(classes);
-
 #   if EXPLICIT_NORMAL_EDGES
     std::vector<bool> edgeList(edges, 1);
 #   endif
 
-    for (auto e : spanningTree)
-      classEdges.push_back(e - 1);
-
     for (const auto& elemCycle : elementaryCycles) {
       nonClassEdges.push_back(std::make_pair(elemCycle[0] - 1, Mapping(edges)));
       Mapping &map = nonClassEdges.back().second;
@@ -108,9 +99,24 @@ struct Tester : public AbstractTester {
 #   endif
 
 #   if EXPLICIT_NORMAL_EDGES
+    for (auto e : spanningTree)
+      if (edgeList[e - 1]) classEdges.push_back(e - 1);
+
     for (size_t i = 0; i < edges; i++)
       if (edgeList[i]) normalEdges.push_back(i);
+#   else
+    for (auto e : spanningTree)
+      classEdges.push_back(e - 1);
+#   endif
+
+    size_t num_classes = 1;
+    edges_ = classEdges.size();
+    while (edges_-- > 0) num_classes *= Ring::size;
+#   if USE_TWO_CUTS
+    for (const auto& cut : twoCuts) num_classes *= cut.num_classes;
 #   endif
+    ClassesType C(num_classes, 0);
+    C.swap(classes);
   }
 
 
@@ -166,23 +172,11 @@ struct Tester : public AbstractTester {
         map.combine(map[i.first], i.second);
     }
 
-#   if USE_TWO_CUTS
-    for (const auto &c : twoCuts) {
-      T m_a = map[c.a];
-      T m_b = map[c.b];
-      if (map[c.a] > map[c.b]) {
-        map.assign(c.a, m_b);
-        map.assign(c.b, m_a);
-      }
-    }
-#   endif
-
     return map;
   }
 
   size_t pack(const Mapping& map) {
     size_t index = 0;
-
     for (auto i : classEdges) {
       index *= Ring::size;
       index += map[i];
@@ -190,8 +184,8 @@ struct Tester : public AbstractTester {
 
 #   if USE_TWO_CUTS
     for (const auto& c : twoCuts) {
-      index *= c.max_value;
-      index += c.getValue(map);
+      index *= c.num_classes;
+      index += c.getClass(map);
     }
 #   endif
 
@@ -223,6 +217,10 @@ struct Tester : public AbstractTester {
 #   if USE_NEXT_FORB
       for (size_t i = 0; i < edges; i++) forb.assign(i, 1);
 
+#     if USE_TWO_CUTS
+      for (const auto& cut : twoCuts) cut.encode(forb, 0);
+#     endif
+
       do {
         Mapping copy(forb);
         INC(classes[pack(canonize(copy))]);

twoCuts.h

@@ -6,25 +6,44 @@ struct TwoCutInt {
   bool flip;
 
   static const int max_value = ((Ring::size - 1) * (Ring::size - 2)) / 2;
+  static const int num_classes = (Ring::size * (Ring::size - 1)) / 2;
 
   struct Tables {
     char decodeMatrix[Ring::size][Ring::size];
     struct { char a, b; } encodeTable[max_value];
+    char classMatrix[Ring::size][Ring::size];
 
     Tables() {
       encodeTable[0] = {1, 2};
 
       for (int i = 1; i < max_value; i++) {
-        bool carry = (encodeTable[i-1].b >= Ring::size);
+        bool carry = (encodeTable[i-1].b + 1 >= Ring::size);
         encodeTable[i].a = encodeTable[i-1].a + carry;
         encodeTable[i].b = carry ? encodeTable[i].a + 1 : encodeTable[i-1].b + 1;
       }
 
       memset(decodeMatrix, -1, sizeof(decodeMatrix));
 
-      for (int i = 0; i < max_value; i++)
+      for (int i = 0; i < max_value; i++) {
         decodeMatrix[encodeTable[i].a][encodeTable[i].b] = i;
       }
+
+      memset(classMatrix, -1, sizeof(classMatrix));
+
+      struct { char a, b; } classTable[num_classes];
+      classTable[0] = {0, 1};
+
+      for (int i = 1; i < num_classes; i++) {
+        bool carry = (classTable[i-1].b + 1 >= Ring::size);
+        classTable[i].a = classTable[i-1].a + carry;
+        classTable[i].b = carry ? classTable[i].a + 1 : classTable[i-1].b + 1;
+      }
+
+      for (int i = 0; i < num_classes; i++) {
+        classMatrix[classTable[i].a][classTable[i].b] = i;
+        classMatrix[classTable[i].b][classTable[i].a] = i;
+      }
+    }
   };
 
   static const Tables S;
@@ -38,6 +57,11 @@ struct TwoCutInt {
     map.assign(a, vals.a);
     map.assign(b, flip ? Ring::negate(vals.b) : vals.b);
   }
+
+  inline int getClass(const Mapping &map) const {
+    int ret = S.classMatrix[map[a]][flip ? Ring::negate(map[b]) : map[b]];
+    return S.classMatrix[map[a]][flip ? Ring::negate(map[b]) : map[b]];
+  }
 };
 
 #endif