minimal change ordering fix

include/permanent/combinatoric.h

@@ -5,18 +5,20 @@
 
 #include <permanent/common.h>
 #include <permanent/complex.h>
-#include <permanent/kperm-gray.h>
-#include <permanent/perm-mv0.h>
+#include <permanent/fxt/kperm-gray.h>
+#include <permanent/fxt/perm-mv0.h>
 
 namespace permanent {
 
 template <typename T, typename I = void>
-result_t<T, I> combinatoric_square(const size_t m, const size_t n, const T *ptr)
+result_t<T, I> combinatoric_square(const size_t m, const size_t, const T *ptr)
 {
-  (void)n;
-
   perm_mv0 permutations(m);
 
+  if (m == 1) {
+    return ptr[0];
+  }
+
   result_t<T, I> out = 0;
   const size_t *perm = permutations.data();
   do {
@@ -58,4 +60,4 @@ result_t<T, I> combinatoric(const size_t m, const size_t n, const T *ptr)
 
 }  // namespace permanent
 
-#endif  // permanent_combinatoric_h_
+#endif  // permanent_combinatoric_h_

include/permanent/common.h

@@ -3,15 +3,17 @@
 #if !defined(permanent_common_h_)
 #define permanent_common_h_
 
-#include <complex>
 #include <cstdint>
+
+#include <complex>
 #include <type_traits>
 
 namespace permanent {
 
 // Numerical types
 
 using size_t = std::size_t;
+using ptrdiff_t = std::ptrdiff_t;
 
 using sgn_t = std::int_fast8_t;
 
@@ -51,4 +53,4 @@ using result_t = _result_t<Type, IntType>::type;
 
 }  // namespace permanent
 
-#endif  // permanent_common_h_
+#endif  // permanent_common_h_

include/permanent/complex.h

@@ -18,21 +18,21 @@ struct _identity_t
 template <typename T>
 using identity_t = _identity_t<T>::type;
 
-}  // namespace permanent
-
-#define DEF_COMPLEX_OPS(OP)                                                                \
-                                                                                           \
-  template <typename T>                                                                    \
-  std::complex<T> operator OP(const std::complex<T> &x, const permanent::identity_t<T> &y) \
-  {                                                                                        \
-    return x OP y;                                                                         \
-  }                                                                                        \
-                                                                                           \
-  template <typename T>                                                                    \
-  std::complex<T> operator OP(const permanent::identity_t<T> &x, const std::complex<T> &y) \
-  {                                                                                        \
-    return x OP y;                                                                         \
-  }
+namespace complex_ops {
+
+#define DEF_COMPLEX_OPS(OP)                                                                     \
+                                                                                                \
+template <typename T>                                                                           \
+inline std::complex<T> operator OP(const std::complex<T> &x, const permanent::identity_t<T> &y) \
+{                                                                                               \
+  return operator OP(x, y);                                                                     \
+}                                                                                               \
+                                                                                                \
+template <typename T>                                                                           \
+inline std::complex<T> operator OP(const permanent::identity_t<T> &x, const std::complex<T> &y) \
+{                                                                                               \
+  return operator OP(x, y);                                                                     \
+}
 
 DEF_COMPLEX_OPS(+)
 DEF_COMPLEX_OPS(-)
@@ -41,4 +41,8 @@ DEF_COMPLEX_OPS(/)
 
 #undef DEF_COMPLEX_OPS
 
-#endif  // PERMANENT_COMPLEX_H_
+} // namespace complex_ops
+
+}  // namespace permanent
+
+#endif  // PERMANENT_COMPLEX_H_

include/permanent/fxt/kperm-gray.h

@@ -0,0 +1,152 @@
+#if !defined HAVE_KPERM_GRAY_H__
+#define HAVE_KPERM_GRAY_H__
+// This file is part of the FXT library.
+// Copyright (C) 2010, 2012, 2013, 2014, 2019, 2023 Joerg Arndt
+// License: GNU General Public License version 3 or later,
+// see the file COPYING.txt in the main directory.
+
+#include <cstdlib>
+
+#include <permanent/fxt/swap.h>
+
+class kperm_gray
+// Gray code for k-permutations of n elements.
+// Same as: k-prefixes of permutations of n elements.
+// Same as: arrangements of k out of n elements.
+// CAT algorithm based on mixed radix Gray code
+//   for the factorial number system (falling base).
+{
+ protected:
+  std::size_t d_[64];      // mixed radix digits with radix = [n-1, n-2, ..., 2]
+  std::size_t i_[64];      // directions
+  std::size_t ix_[64];     // permutation (inverse perms in Trotter's order)
+  std::size_t x_[64];      // inverse permutation (perms in Trotter's order)
+  std::size_t n_;          // total of n elements
+  std::size_t k_;          // prefix length: permutations of k elements
+  std::size_t sw1_, sw2_;  // indices of elements swapped most recently
+
+  kperm_gray(const kperm_gray &) = delete;
+  kperm_gray &operator=(const kperm_gray &) = delete;
+
+ public:
+  explicit kperm_gray(std::size_t n)
+  {
+    n_ = n;
+    k_ = n;
+    // d_ = new std::size_t[n_];
+    d_[n - 1] = -1UL;  // sentinel
+    // i_ = new std::size_t[n_];
+    // x_ = new std::size_t[n_];
+    // ix_ = new std::size_t[n_];
+    i_[n_ - 1] = 0;
+    first(n_);
+  }
+
+  ~kperm_gray()
+  {
+    // delete [] i_;
+    // delete [] d_;
+    // delete [] x_;
+    // delete [] ix_;
+  }
+
+  const std::size_t *data() const { return ix_; }
+  const std::size_t *invdata() const { return x_; }
+  void get_swap(std::size_t &s1, std::size_t &s2) const
+  {
+    s1 = sw1_;
+    s2 = sw2_;
+  }
+
+  const std::size_t *mr_digits() const { return d_; }
+
+  void first(std::size_t k)
+  {
+    k_ = k;
+    for (std::size_t j = 0; j < n_ - 1; ++j) d_[j] = 0;
+    for (std::size_t j = 0; j < n_ - 1; ++j) i_[j] = +1;
+    for (std::size_t j = 0; j < n_; ++j) x_[j] = ix_[j] = j;
+    sw1_ = n_ - 1;
+    sw2_ = n_ - 2;
+  }
+
+  //    void last(std::size_t k)
+  //    {
+  //        first(k);  // k_ is set with call
+  //        d_[n_-2] = 1;
+  //        swap2(x_[n_-1], x_[n_-2]);
+  //        swap2(ix_[n_-1], ix_[n_-2]);
+  //        for (std::size_t j=0; j<n_-2; ++j)  i_[j] = -1UL;
+  //        i_[n_-2] = +1;
+  //    }
+
+ private:
+  void swap(std::size_t j, std::size_t im)
+  {
+    const std::size_t x1 = j;        // element j
+    const std::size_t i1 = ix_[x1];  // position of j
+    const std::size_t i2 = i1 + im;  // neighbor
+    const std::size_t x2 = x_[i2];   // position of neighbor
+    x_[i1] = x2;
+    x_[i2] = x1;  // swap2(x_[i1], x_[i2]);
+    ix_[x1] = i2;
+    ix_[x2] = i1;  // swap2(ix_[x1], ix_[x2]);
+    sw1_ = i1;
+    sw2_ = i2;
+  }
+
+ public:
+  bool next()
+  {
+    std::size_t j = 0;
+    std::size_t m1 = n_ - 1;  // nine in falling factorial base
+    std::size_t ij;
+    while ((ij = i_[j])) {
+      std::size_t im = i_[j];
+      std::size_t dj = d_[j] + im;
+      if (dj > m1)  // =^= if ( (dj>m1) || ((long)dj<0) )
+      {
+        i_[j] = -ij;
+      } else {
+        d_[j] = dj;
+        swap(j, im);
+        return true;
+      }
+
+      --m1;
+      ++j;
+      if (j >= k_) return false;
+    }
+    return false;
+  }
+
+  //    bool prev()
+  //    {
+  //        std::size_t j = 0;
+  //        std::size_t m1 = n_ - 1;
+  //        std::size_t ij;
+  //        while ( (ij=i_[j]) )
+  //        {
+  //            std::size_t im = -i_[j];
+  //            std::size_t dj = d_[j] + im;
+  //            if ( dj>m1 )  // =^= if ( (dj>m1) || ((long)dj<0) )
+  //            {
+  //                i_[j] = -ij;
+  //            }
+  //            else
+  //            {
+  //                d_[j] = dj;
+  //                swap(j, im);
+  //                return true;
+  //            }
+  //
+  //            --m1;
+  //            ++j;
+  //            if ( j>=k_ )  return false;
+  //        }
+  //        return false;
+  //    }
+};
+// -------------------------
+
+#endif  // !defined HAVE_KPERM_GRAY_H__

include/permanent/fxt/perm-mv0.h

@@ -0,0 +1,120 @@
+#if !defined HAVE_PERM_MV0_H__
+#define HAVE_PERM_MV0_H__
+// This file is part of the FXT library.
+// Copyright (C) 2010, 2011, 2012, 2014, 2019, 2023 Joerg Arndt
+// License: GNU General Public License version 3 or later,
+// see the file COPYING.txt in the main directory.
+
+#include <cstdlib>
+
+#include <permanent/fxt/swap.h>
+
+// define to update d[0] with each step:
+#define PERM_MV0_UPDATE_D0  // default is on
+
+// whether to use arrays instead of pointers:
+#define PERM_MV0_FIXARRAYS  // small speedup; default off
+
+class perm_mv0
+// Inverse permutations corresponding to falling factorial numbers.
+// CAT algorithm based on mixed radix Gray code
+//   for the factorial number system (falling base).
+{
+ public:
+#if !defined PERM_MV0_FIXARRAYS
+  std::size_t *d_;  // mixed radix digits with radix = [n-1, n-2, n-3, ..., 2]
+  std::size_t *x_;  // permutation
+#else
+  std::size_t d_[64];
+  std::size_t x_[64];
+#endif
+  std::size_t ect_;  // counter for easy case
+  std::size_t n_;    // permutations of n elements
+
+  perm_mv0(const perm_mv0 &) = delete;
+  perm_mv0 &operator=(const perm_mv0 &) = delete;
+
+ public:
+  explicit perm_mv0(std::size_t n)
+  // Must have n>=2
+  {
+    n_ = n;
+#if !defined PERM_MV0_FIXARRAYS
+    d_ = new std::size_t[n_];
+    x_ = new std::size_t[n_];
+#endif
+    d_[n - 1] = 1;  // sentinel (must be nonzero)
+    first();
+  }
+
+  ~perm_mv0()
+  {
+#if !defined PERM_MV0_FIXARRAYS
+    delete[] d_;
+    delete[] x_;
+#endif
+  }
+
+  const std::size_t *data() const { return x_; }
+
+  void first()
+  {
+    for (std::size_t k = 0; k < n_; ++k) x_[k] = k;
+    for (std::size_t k = 0; k < n_ - 1; ++k) d_[k] = 0;
+    ect_ = 0;
+  }
+
+  bool next()
+  {
+    if (++ect_ < n_)  // easy case
+    {
+#if defined PERM_MV0_UPDATE_D0
+      d_[0] = ect_;
+#endif
+      swap2(x_[ect_], x_[ect_ - 1]);
+      return true;
+    } else {
+      ect_ = 0;
+#if defined PERM_MV0_UPDATE_D0
+      d_[0] = ect_;
+#endif
+
+      std::size_t j = 1;
+      std::size_t m1 = n_ - 2;  // nine in falling factorial base
+      while (d_[j] == m1)       // find digit to increment
+      {
+        d_[j] = 0;
+        --m1;
+        ++j;
+      }
+
+      if (j == n_ - 1) return false;  // current permutation is last
+
+      const std::size_t dj = d_[j];
+      d_[j] = dj + 1;
+
+      // element at d[j] moves one position to the right:
+      swap2(x_[dj], x_[dj + 1]);
+
+      {  // move n-j elements to end:
+        std::size_t s = n_ - j, d = n_;
+        do {
+          --s;
+          --d;
+          x_[d] = x_[s];
+        } while (s);
+      }
+
+      // fill in 0,1,2,..,j-1 at start:
+      for (std::size_t k = 0; k < j; ++k) x_[k] = k;
+
+      return true;
+    }
+  }
+};
+// -------------------------
+
+// #undef PERM_MV0_UPDATE_D0
+// #undef PERM_MV0_FIXARRAYS
+
+#endif  // !defined HAVE_PERM_MV0_H__