Fixed Ryser, updated tuning to 8 parameters and two hyperplanes (#25)

* fixed both ryser algs, added tests, updated tuning
* updated variable init, used popcount, default tuning off

Makefile

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

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,30 @@
 #include <permanent/tables.h>
 
 #include <bit>
+#include <vector>
 
 namespace {
 
+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
+  k_perms.reserve(1 << n);
+
+  for (int i = 0; i < (1 << n); i++) {
+    unsigned int cur = i ^ (i >> 1);
+    if (std::popcount(cur) == k) {
+      std::vector<int> vec;  // Initialize vector for new combination
+      vec.reserve(n);
+      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)
 {
@@ -93,122 +114,65 @@ bool gen_next_perm(T *const falling_fact, T *const perm, T *const inv_perm, cons
 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)
+result_t<T, I> ryser_square(const size_t m, const size_t, 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.
-   */
-
+  T out = 0;
   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.  */
-
+  for (size_t k = 0; k < c; ++k) {
+    T rowsumprod = 1.0;
+    for (size_t i = 0; i < m; ++i) {
+      T rowsum = 0.0;
+      for (size_t j = 0; j < m; ++j) {
         if (k & (1UL << j)) {
           rowsum += ptr[m * i + j];
         }
       }
-
-      /* Update product of row sums. */
-
       rowsumprod *= rowsum;
     }
-
-    /* Add term multiplied by the parity of the characteristic vector. */
-
-    out += rowsumprod * parity(std::popcount(k));
+    out += rowsumprod * (1 - ((std::popcount(k) & 1) << 1));
   }
 
-  /* 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);
-
-    /* sort up to position u + 1 where u = min(m - k, n - 1). */
-    size_t sort_up_to = n - 1;
-
-    if ((m - k) < sort_up_to) {
-      sort_up_to = (m - k);
-    }
-
-    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]];
+  int sign = 1;
+  result_t<T, I> out = 0.0;
+  std::vector<std::vector<int>> k_perms;
+
+  /* Iterate over subsets from size 0 to size m */
+
+  for (size_t k = 0; k < m; ++k) {
+    all_combinations(n, m - k, k_perms);
+    size_t bin = binomial((n - m + k), k);
+    result_t<T, I> permsum = 0.0;
+
+    /* Compute permanents of each submatrix */
+    result_t<T, I> vec[64];
+    for (const auto &combination : k_perms) {
+      result_t<T, I> colprod;
+      for (size_t i = 0; i < m; ++i) {
+        result_t<T, I> matsum = 0.0;
+        for (size_t j = 0; j < (m - k); ++j) matsum += ptr[i * n + combination[j]];
+        vec[i] = matsum;
       }
-      vec[i] = sum_of_matrix_vals;
-    }
-    for (i = 0; i < m; ++i) {
-      prod_of_cols *= vec[i];
+
+      colprod = 1.0;
+      for (size_t i = 0; i < m; ++i) colprod *= vec[i];
+      permsum += colprod * sign * bin;
     }
 
-    result += value_sign * bin_c * prod_of_cols;
-    counter += 1;
+    /* Add term to result */
 
-    /* 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];
-      }
+    out += permsum;
 
-      result += value_sign * bin_c * prod_of_cols;
-      counter += 1;
-    }
-    sum_over_k_vals += (result / counter) * binomial(n, (m - k));
-    value_sign *= -1;
+    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

@@ -37,6 +37,7 @@ test = ["pytest"]
 
 [tool.scikit-build]
 wheel.expand-macos-universal-tags = true
+cmake.args = ["-DPERMANENT_TUNE=OFF"]
 
 [tool.pytest.ini_options]
 minversion = "6.0"