fixed both ryser algs, added tests, updated tuning

Makefile

@@ -26,9 +26,9 @@ clean:
 install: _build
 	@cmake --install build
 
+build:
+	@cmake -B build $(CMAKE_FLAGS)
+
 _build: build
 	@cmake --build build --target all
-	@test -e build/compile_commands.json && ln -sf build/compile_commands.json .
-
-build:
-	@cmake -B build $(CMAKE_FLAGS) .
+	@if test -e build/compile_commands.json; then ln -sf build/compile_commands.json .; fi

include/permanent/opt.h

@@ -20,22 +20,49 @@ namespace permanent {
 template <typename T, typename I = void>
 result_t<T, I> opt_square(const size_t m, const size_t n, const T *ptr)
 {
-  if (n <= PARAM_4<T, I>) {
-    return ryser_square<T, I>(m, n, ptr);
-  } else if (n * PARAM_1<T, I> + PARAM_2<T, I> + PARAM_3<T, I> > 0) {
-    return combinatoric_square<T, I>(m, n, ptr);
-  } else {
-    return glynn_square<T, I>(m, n, ptr);
+  if (n <= PARAM_8<T, I>) {  // N ≤ 13
+    if (n <= PARAM_4<T, I>) {
+      return combinatoric_square<T, I>(m, n, ptr);
+    } else {
+      // Check first hyperplane: -6.67(m/n) + (-0.35)(n) + 6.49 = 0
+      // Using stored parameters for flexibility
+      const double ratio = static_cast<double>(m) / n;
+      if (PARAM_1<T, I> * ratio + PARAM_2<T, I> * n + PARAM_3<T, I> > 0) {
+        return combinatoric_square<T, I>(m, n, ptr);
+      } else {
+        return glynn_square<T, I>(m, n, ptr);
+      }
+    }
+  } else {  // N > 13
+    // Check second hyperplane: -10.76(m/n) + 0.09(n) + 1.96 = 0
+    const double ratio = static_cast<double>(m) / n;
+    if (PARAM_5<T, I> * ratio + PARAM_6<T, I> * n + PARAM_7<T, I> > 0) {
+      return glynn_square<T, I>(m, n, ptr);
+    } else {
+      return ryser_square<T, I>(m, n, ptr);
+    }
   }
 }
 
 template <typename T, typename I = void>
 result_t<T, I> opt_rectangular(const size_t m, const size_t n, const T *ptr)
 {
-  if (n * PARAM_1<T, I> + PARAM_2<T, I> * m / n + PARAM_3<T, I> > 0) {
-    return combinatoric_rectangular<T, I>(m, n, ptr);
-  } else {
-    return glynn_rectangular<T, I>(m, n, ptr);
+  if (n <= PARAM_8<T, I>) {  // N ≤ 13
+    // First hyperplane for small matrices
+    const double ratio = static_cast<double>(m) / n;
+    if (PARAM_1<T, I> * ratio + PARAM_2<T, I> * n + PARAM_3<T, I> > 0) {
+      return combinatoric_rectangular<T, I>(m, n, ptr);
+    } else {
+      return glynn_rectangular<T, I>(m, n, ptr);
+    }
+  } else {  // N > 13
+    // Second hyperplane for large matrices
+    const double ratio = static_cast<double>(m) / n;
+    if (PARAM_5<T, I> * ratio + PARAM_6<T, I> * n + PARAM_7<T, I> > 0) {
+      return glynn_rectangular<T, I>(m, n, ptr);
+    } else {
+      return ryser_rectangular<T, I>(m, n, ptr);
+    }
   }
 }
 
@@ -47,4 +74,4 @@ result_t<T, I> opt(const size_t m, const size_t n, const T *ptr)
 
 }  // namespace permanent
 
-#endif  // permanent_opt_h_
+#endif  // permanent_opt_h_

include/permanent/ryser.h

@@ -8,9 +8,40 @@
 #include <permanent/tables.h>
 
 #include <bit>
+#include <vector>
 
 namespace {
 
+/* Helper functions for rectangle ryser. */
+inline int gray_code (int n) {
+    return n ^ (n >> 1);
+}
+
+inline int count_bits (int n) {
+    int res = 0;
+    for (; n; n >>= 1)
+        res += n & 1;
+    return res;
+}
+
+inline void all_combinations(int n, int k, std::vector<std::vector<int> >& k_perms) {
+    k_perms.clear(); // Clear k_perms to make sure it starts empty
+
+    for (int i = 0; i < (1 << n); i++) {
+        int cur = gray_code(i);
+        if (count_bits(cur) == k) {
+            std::vector<int> vec; // Initialize vector for new combination
+
+            for (int j = 0; j < n; j++) {
+                if (cur & (1 << j)) {
+                    vec.push_back(j); // Add element to vector
+                }
+            }
+            k_perms.push_back(vec); // Store the combination
+        }
+    }
+}
+
 template <typename T>
 auto parity(const T &x)
 {
@@ -95,120 +126,72 @@ namespace permanent {
 template <typename T, typename I = void>
 result_t<T, I> ryser_square(const size_t m, const size_t n, const T *ptr)
 {
-  (void)n;
-
-  size_t i, j;
-  size_t k;
-  result_t<T, I> rowsum;
-  result_t<T, I> out = 0;
-
-  /* Iterate over c = pow(2, m) submatrices (equal to (1 << m)) submatrices.
-   */
-
-  size_t c = 1UL << m;
-
-  /* Loop over columns of submatrix; compute product of row sums. */
-
-  for (k = 0; k < c; ++k) {
-    result_t<T, I> rowsumprod = 1.0;
-    for (i = 0; i < m; ++i) {
-      /* Loop over rows of submatrix; compute row sum. */
-
-      rowsum = 0.0;
-      for (j = 0; j < m; ++j) {
-        /* Add element to row sum if the row index is in the
-         * characteristic * vector of the submatrix, which is the binary
-         * vector given by k.  */
-
-        if (k & (1UL << j)) {
-          rowsum += ptr[m * i + j];
+    (void)n;
+
+    size_t i, j, k;
+    T rowsum;  // Use T instead of result_t for intermediate calculations
+    T out = 0;
+    size_t c = 1UL << m;
+
+    for (k = 0; k < c; ++k) {
+        T rowsumprod = 1.0;
+        for (i = 0; i < m; ++i) {
+            rowsum = 0.0;
+            for (j = 0; j < m; ++j) {
+                if (k & (1UL << j)) {
+                    rowsum += ptr[m * i + j];
+                }
+            }
+            rowsumprod *= rowsum;
         }
-      }
-
-      /* Update product of row sums. */
-
-      rowsumprod *= rowsum;
+        out += rowsumprod * (1 - ((__builtin_popcountll(k) & 1) << 1));
     }
 
-    /* Add term multiplied by the parity of the characteristic vector. */
-
-    out += rowsumprod * parity(std::popcount(k));
-  }
-
-  /* Return answer with the correct sign (times -1 for odd m). */
-
-  return out * parity(m);
+    return static_cast<result_t<T, I>>(out * ((m % 2 == 1) ? -1 : 1));
 }
 
 template <typename T, typename I = void>
 result_t<T, I> ryser_rectangular(const size_t m, const size_t n, const T *ptr)
 {
-  size_t falling_fact[128];
-  size_t perm[128];
-  size_t inv_perm[128];
-  falling_fact[0] = 0;
-
-  size_t i, j, k;
-
-  int value_sign = 1;
-
-  result_t<T, I> sum_of_matrix_vals;
-  result_t<T, I> sum_over_k_vals = 0.0;
-  result_t<T, I> vec[64];
-
-  /* Dealing with a rectangle. Can't use bit hacking trick here. */
-  for (k = 0; k < m; ++k) {
-    /* Store the binomial coefficient for this k value bin_c. */
-    size_t counter = 0;  // Count how many permutations you have generated
-    size_t bin_c = binomial((n - m + k), k);
-    result_t<T, I> prod_of_cols = 1;
-    result_t<T, I> result = 0;
-
-    /* (Re)initialize the set to permute for this k value. */
-    init_perm(n, falling_fact, perm, inv_perm);
+    size_t i, j, k, bin;
+    int sign = 1;
+
+    result_t<T, I> colprod, matsum, permsum;
+    result_t<T, I> out = 0.0;
+    result_t<T, I> vec[64];
+    std::vector<std::vector<int> > k_perms;
+
+    /* Iterate over subsets from size 0 to size m */
+
+    for (k = 0; k < m; ++k) {
+        all_combinations(n, m - k, k_perms);
+
+        bin = binomial((n - m + k), k); 
+        permsum = 0.0;
+
+        /* Compute permanents of each submatrix */
+        for (const auto& combination : k_perms) {
+                for (i = 0; i < m; ++i) {
+                    matsum = 0.0;
+                    for (j = 0; j < (m - k); ++j)
+                        matsum += ptr[i * n + combination[j]];
+                    vec[i] = matsum;
+                }
+
+                colprod = 1.0;
+                for (i = 0; i < m; ++i)
+                    colprod *= vec[i];
+                permsum += colprod * sign * bin;
+        }
 
-    /* sort up to position u + 1 where u = min(m - k, n - 1). */
-    size_t sort_up_to = n - 1;
+        /* Add term to result */
 
-    if ((m - k) < sort_up_to) {
-      sort_up_to = (m - k);
-    }
+        out += permsum;
 
-    for (i = 0; i < m; ++i) {
-      sum_of_matrix_vals = 0.0;
-      for (j = 0; j < sort_up_to; ++j) {
-        sum_of_matrix_vals += ptr[i * n + perm[j]];
-      }
-      vec[i] = sum_of_matrix_vals;
-    }
-    for (i = 0; i < m; ++i) {
-      prod_of_cols *= vec[i];
+        sign *= -1;
     }
 
-    result += value_sign * bin_c * prod_of_cols;
-    counter += 1;
-
-    /* Iterate over second to last permutations of the set. */
-    while (gen_next_perm(falling_fact, perm, inv_perm, n, sort_up_to)) {
-      for (i = 0; i < m; ++i) {
-        sum_of_matrix_vals = 0.0;
-        for (j = 0; j < m - k; ++j) {
-          sum_of_matrix_vals += ptr[i * n + perm[j]];
-        }
-        vec[i] = sum_of_matrix_vals;
-      }
-      prod_of_cols = 1.0;
-      for (i = 0; i < m; ++i) {
-        prod_of_cols *= vec[i];
-      }
-
-      result += value_sign * bin_c * prod_of_cols;
-      counter += 1;
-    }
-    sum_over_k_vals += (result / counter) * binomial(n, (m - k));
-    value_sign *= -1;
-  }
-  return sum_over_k_vals;
+    return out;
 }
 
 template <typename T, typename I = void>

include/permanent/tuning.default.h

@@ -8,10 +8,17 @@ namespace permanent {
 template <typename Type, typename IntType = void>
 struct _tuning_params_t
 {
-  static constexpr double PARAM_1 = -5.72098e-01;
-  static constexpr double PARAM_2 = -2.20142e+01;
-  static constexpr double PARAM_3 = +1.52979e+01;
-  static constexpr double PARAM_4 = +3.00000e+00;
+  // First hyperplane parameters (N ≤ 13): -6.67(m/n) + (-0.35)(n) + 6.49 = 0
+  static constexpr double PARAM_1 = -6.67000e+00;  // coefficient of m/n
+  static constexpr double PARAM_2 = -3.50000e-01;  // coefficient of n
+  static constexpr double PARAM_3 = +6.49000e+00;  // constant term
+  static constexpr double PARAM_4 = +4.00000e+00;  // Combinatorial crossover limit
+
+  // Second hyperplane parameters (N > 13): -10.76(m/n) + 0.09(n) + 1.96 = 0
+  static constexpr double PARAM_5 = -1.07600e+01;  // coefficient of m/n
+  static constexpr double PARAM_6 = +9.00000e-02;  // coefficient of n
+  static constexpr double PARAM_7 = +1.96000e+00;  // constant term
+  static constexpr double PARAM_8 = +1.30000e+01;  // Combinatorial limit (N=13)
 };
 
 template <typename Type, typename IntType = void>
@@ -26,6 +33,18 @@ static constexpr double PARAM_3 = _tuning_params_t<Type, IntType>::PARAM_3;
 template <typename Type, typename IntType = void>
 static constexpr double PARAM_4 = _tuning_params_t<Type, IntType>::PARAM_4;
 
+template <typename Type, typename IntType = void>
+static constexpr double PARAM_5 = _tuning_params_t<Type, IntType>::PARAM_5;
+
+template <typename Type, typename IntType = void>
+static constexpr double PARAM_6 = _tuning_params_t<Type, IntType>::PARAM_6;
+
+template <typename Type, typename IntType = void>
+static constexpr double PARAM_7 = _tuning_params_t<Type, IntType>::PARAM_7;
+
+template <typename Type, typename IntType = void>
+static constexpr double PARAM_8 = _tuning_params_t<Type, IntType>::PARAM_8;
+
 }  // namespace permanent
 
-#endif  // permanent_tuning_h_
+#endif  // permanent_tuning_h_

pyproject.toml

@@ -38,6 +38,9 @@ test = ["pytest"]
 [tool.scikit-build]
 wheel.expand-macos-universal-tags = true
 
+[tool.scikit-build]
+cmake.args = ["-DPERMANENT_TUNE=ON"]
+
 [tool.pytest.ini_options]
 minversion = "6.0"
 addopts = ["-ra", "--showlocals", "--strict-markers", "--strict-config"]